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New Thermomechanical Reciprocity Relations 
With Application to Thermal Stress Analysis’ 


M. A. BIOT* 


Cornell Aeronautical Laboratory, Inc. 


SUMMARY 


Based on the variational formulation of linear thermodynam- 
ics as developed previously by the writer, thermomechanical 
reciprocity relations are discussed which lead to new methods of 
analysis of thermal stresses. These reciprocity relations are quite 
different from the usual ones derived from the analogy of thermal 
loading with a combination of surface and body-force distribution, 
The results are applicable to stationary and transient temper- 
atures in elastic and viscoelastic structures. The methods are 
entirely variational and do not require the evaluation of the 
temperature field. The stresses at one point are expressed di- 
rectly in terms of any arbitrary distribution temperatures applied 
externally, including the effect cf surface heat-transfer layer. 
The concepts and procedures are illustrated on a simple ex- 
ample. The relation is pointed out between the reciprocity 
property and the generalization of Castigliano’s principle to 
thermomechanics. 


(1) INTRODUCTION 


N THE PAST, problems of thermal stresses in elastic 
I systems have been treated by neglecting the 
reciprocal coupling between the temperature and 
deformation. What is meant by reciprocal coupling 
is the fact that a change of temperature produces 
a deformation and in turn a deformation produces 
a change of temperature. Classical thermodynamics 
shows that one effect cannot occur without the other. 
This coupling leads to the well-known phenomenon of 
thermoelastic dissipation in elastic solids. The cou- 
pling is important from the standpoint of the physicist 
and in certain specific technological applications such 
as electronics, but in the theory of structures its order 
of magnitude is not significant. 

However the coupling does acquire importance in 
thermal stress analysis as a concept because, as we will 
show, it leads to entirely new methods of calculation. 
We have shown? that the most general thermoelastic 
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system obeys variational principles. These variational 
principles were derived by introducing the concept of 
entropy displacement and thermoelastic potential and 
by expressing a dissipation function in terms of the time 
derivatives of the entropy displacement. From these 
principles we have derived a formulation of thermo- 
elasticity where the mechanical and thermal variables 
play identical roles. The boundary forces and tem- 
peratures then obey thermomechanical reciprocity rela- 
tions. In the case of stationary temperatures, this is 
analogous to Maxwell’s reciprocity relations for the 
forces acting on an elastic structure. Furthermore, 
the reciprocity properties apply to transient temper- 
atures. Attention is called to the essential difference 
between the reciprocity relations discussed in this 
paper and the usual ones derived from the analogy of 
thermal loading with distributed surface and body 
forces or its equivalent variational form which uses the 
isothermal free energy. The new viewpoint presented 
here leads to influence coefficients and influence func- 
tions which are truly of a hybrid thermomechanical char- 
acter -i.e., they belong to both thermodynamics and 
mechanics. By their use it becomes possible to predict 
the stresses of one point due to any application of tem- 
perature at the boundary applied either directly to the 
solid or through a surface heat-transfer layer. Ad- 
vantages of the method are many. The procedure 
makes use entirely of variational methods and by- 
passes the evaluation of the temperature field itself. The 
process is error-smoothing, eliminates the evaluation of 
eventual local temperature singularities, and requires 
only one calculation for all possible cases of applied 
temperatures at the boundary. The procedure may 
also be used to calculate thermal stresses in visco- 
elastic systems. 

In Section (2) we review our previous results on the 
variational formulation of thermoelasticity. The re- 
sulting reciprocity relations for stationary temper- 
atures are discussed in Section (3), along with some 
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variational procedures for the evaluation of the thermo- 
mechanical influence coefficients. Section (4) dis- 
cusses the use of the reciprocity relations in thermal 
stress analysis for stationary temperatures, and Sec- 
tion (5) extends the same methods for transient thermal 
stresses in elastic and viscoelastic structures. For the 
purpose of illustration of the concepts and procedure, 
we have treated a simple example in Section (6). 


(2) VARIATIONAL FORMULATION OF THERMOELASTICITY 


In previous work*?~® we have described the state of a 
thermoelastic system by two vector fields u and S of 
components u; and S,, respectively. The vector u is 
the geometrical displacement of the medium. The 
vector S, which we have called the entropy displace- 
ment or entropy flow, is defined so that its time deriv- 
ative is 


(0S/ot) = (1/T,) (OH/Ot) (2.1) 


where 0H/Ot is the rate of heat flow. The symbol 7, 
denotes a constant reference temperature which repre- 
sents the uniform temperature acquired by the system 
when it is in a state of equilibrium in the presence of a 
large heat reservoir at the same temperature 7,. The 
terminology “entropy displacement” for S is of course 
only justified if the local temperatures 6 + 7, are such 
that@<<T7;,. The theory, being linear, is valid when 
this inequality is not fulfilled, but there is no need of 
changing the terminology. The thermoelastic proper- 
ties of the system are then completely defined’ by a 
thermoelastic potential 


V = ff [W + (1/2) (T,/c) (Bye + div S)*]dr 
(2.2) 


and a dissipation function 


D = (1/2)T, fff (0.S;/Ot)dr + 


(1/2)T, ff. (1/K) (0S,/0t)*dA (2.3) 


The integrals are extended to the volume 7 and the 
boundary surface A. The symbols are 


ey = (1/2) [(Ou;/dx;) + (Ou;/Ox;)] strain com- 
ponents 
C,.'’ = isothermal elastic moduli 
W isothermal strain energy 
K = coefficient of surface heat transfer 
= = thermal resistivity matrix—i.e., 
inverse of the thermal conductivity (k;,) 
‘ S, = normal component of S at the boundary 
c = specific heat per unit volume at zero strain 
a The §;;’s are coefficients which appear in the equations 
of state for the stress 
Cuy = By (2.4) 


where @ is the local excess temperature above the 
equilibrium value 7,. The field S, the strain e,,, and 
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the excess temperature @ are related by the equation 
div Ss = —(c6/T,) ij (2.5) 


We now disturb the equilibrium of the system by 

applying forces F per unit area at the boundary and 
external temperatures 0,. These temperatures may be 
applied at the boundary outside the surface heat-transfer 
layer and not directly to the solid. They may also of 
course be applied directly to the solid, which amounts 
to putting K = o. 
We have shown that the system responds as if general- 
ized thermomechanical forces Q; were applied to it, and 
that these forces are defined by a method of virtual 
work. This virtual work is 


= ff (F-du + 0,6S,)dA (2.6) 
A 


where 5S, is the normal component of S at the boundary 
chosen positive inward. The thermoelastic system is 
represented by the generalized coordinates q;, i.e., we 
have put 

where u; and S, are vector fields of fixed configuration. 
The thermoelastic potential and the dissipation func- 
tion then become 


D = (1/2)bigigs (2.8) 


Applying variational methods, we have derived the 
following differential equations: 


(OV/0q;) + = Q: (2.9) 
or + bigs = Q: 


V= 


These equations are a particular case of the more gen- 
eral treatment of linear thermodynamics which we have 
introduced in references 1 and 2 and displayed in more 
detail in reference 5. The quadratic forms V and D 
are positive definite because of their physical nature. 


(3) RECIPROCITY RELATIONS FOR THE THERMOELASTIC 
ADMITTANCE 


It is of interest to consider the system as a ‘‘black 
box,’’ forces being applied to a small number of “‘ob- 
served”’ coordinates g;. The other coordinates which 
may be very large or even infinite in number are un- 
observed and the conjugate forces applied to them 
are zero. We may solve the system (2.9) for the 
observed coordinates in terms of the corresponding 
forces. We have shown that the solution is 


= Ay*Q, (3.1) 


where A;,;* = A;,;* is a symmetric matrix representing 
the thermoelastic admittance of the system. In oper- 
ational form this admittance is 


= +A] + Cy (3.2) 


This was given a rigorous proof in reference 1. The 
symbol p refers to the time derivative p = d/dt which 


in 
quel 
cour 
tran 
and 

tem] 


prod 


We 
(3.1) 
titie: 
the t 


prod 
whicl 


Asa 
matr. 


This 
mom 
relati 
locati 
plays 
that | 
surfac 
For 
coeffic 
value 
Since 
defor: 
Howe 
eleme 
tain a 
which 


In thi 
tion ¢ 
may | 
slow 1 
the ra 
ciples 
minin 


< 
slow 
bein; 
| with 
Wer 
be z 
finite 
We. 
a area he 


IC 


NEW THERMOMECHANICAL RECIPROCITY RELATIONS 403 


in case of harmonic functions of time with angular fre- 
quency w becomes p = ww. Relations (3.1) may of 
course also be considered as relations between Laplace 
transforms of g; and Q;. Consider now two points P 
and M at the boundary of a thermoelastic system. A 
temperature 6p applied externally at P over a unit area 


produces at M a displacement uy in a certain direction. t 


We may assimilate the subscripts 7 and 7 in expressions 
(3.1) and (3.2) with the points P and M and the quan- 
tities Q; and g; with 6p and uy. We shall assume that 
the temperatures are varying with time at an infinitely 
slow rate, so that at every moment the temperatures 
are the same as in a steady state. This assumption 
being equivalent to that of zero frequency we may put 
p = 0 in expression (3.2). Hence, we may write 


uy = (3.3) 


with an “‘influence coefficient”’ 


Amp = (Cy/ds) + Cy (3.4) 


We note that, if Cc" ~ 0, the corresponding A, cannot 
be zero, since a finite temperature must produce a 
finite displacement. Conversely, a force Fy applied 
at M in the direction uy at an infinitely slow rate will 
produce per unit area at P an inflow of entropy Sp, 
which may be written 


Sp = Apu Fu (3.5) 


As a consequence of the symmetry of the admittance 
matrix (3.2), we may write 


= (3.6) 


This reciprocity relation is the extension to mixed ther- 
momechanical variables of the analogous Maxwell’s 
relations between forces and displacements at different 
locations of an elastic system. The temperature @p 
plays the role of 0,, and we have already pointed out 
that this indicates the temperature either outside the 
surface transfer layer or at the solid boundary itself. 
For practical purposes an interesting feature of the 
coefficient Apy is the possibility of establishing its 
value by a relatively simple procedure, as follows. 
Since the force is applied gradually and very slowly the 
deformation may be considered isothermal (@ = 0). 
However, a flow of entropy is produced because each 
element of volume due to its deformation exudes a cer- 
tain amount of heat. This can be seen from Eq. (2.5), 
which for isothermal deformations, i.e., 8 = 0, becomes 


div S = —Biyli (3.7) 


In this equation e;; is the strain produced by the applica- 
tion of the force Fy. The force Fy and the vector S 
may be assumed to vary linearly with time at a very 
slow rate. The time derivative of S is proportional to 
the rate of heat flow and, according to the general prin- 
ciples developed earlier, must satisfy a principle of 
minimum dissipation—i.e., its spatial distribution must 


+ Note that uy has the dimension of a displacement per unit 
area hence of (length)=}. 


minimize the dissipation function as given by Eq. 
(2.3). Since we are dealing with a stationary state, 
S may be assumed to vary linearly with time and we 
may replace 0S/Ot by S in the dissipation function D. 
The problem is then to minimize 


= (1/2) fff dij + 
(1/2) ff (1/K) (S")*dA (3.8) 
A 


under the constraint (3.7). To do this we may either 


minimize absolutely 
+ ff (div Ss + By€y)dr (3.9) 


with the Lagrangian multiplier \, or, preferably, intro- 
duce any field S* which satisfies Eq. (3.7), then put 


S = S++ S* (3.10) 


where S* is an unknown field with conservative flow 
i.€., 


div S+ = 0 (3.11) 


We then find the field S*+ which minimizes D’ abso- 
lutely. Identical procedures have been discussed by 
this writer in more detail in connection with heat flow 
analysis problems.‘ It may be convenient to express 
the right-hand side of Eq. (3.7) in terms of the stresses 
due to Fy. This will be the case in particular if such 
stresses are statically determined directly from Fy. 
Consider a material under zero stress ¢,, = 0. Rela- 
tion (2.4) becomes 


C,,“ey = (3.12) 


With thermal dilation coefficient a,; we also have 
Cy a9 (3.13) 
C,.%ay = By» (3.14) 


If now e,, is the strain due to the force Fy for isothermal 
deformation (6 = 0), the corresponding stress is 


hence, 


oy = Chey (3.15) 
Multiplying Eq. (3.14) by e,, and taking Eq. (3.15) 
into account we derive 


= pve py (3. 16) 
We may therefore write Eq. (3.7) as 
div = (3.17) 


In the particular case of an isotropic medium, the stress- 
strain relations are 


Cy = Que ij + (Ae B0)5;; 


where € = @;; + Cy, + €z is the dilation, and yu, A, the 
Lamé constants. The coefficients 8;; are reduced to 


Bi = (3.19) 


and 8 may be expressed in terms of the linear coefficient 
a, of thermal dilation as 


(3.18) 
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Fic. 1. Thermal stress at J due to temperature 6P applied at P. 


B = a(2u + 3d) (3.20) 


Relation (3.7) becomes 


div S = —6e (3.21) 


where e¢ is the dilation produced by the force Fy,. This 
equation may also be written in terms of the stress in- 
variant o = o;; + o,, + o., due to Fy by applying 
Eq. (3.17). We find 


div S = —a,a (3.22) 


This expression will be convenient in problems where 
the stress field due to Fy, is statically determined. 
The expression for D’ also simplifies to 


D’ = (1/2) ff fc + 
(1/2) (1/K)S,°dA (3.23) 
A 


where & is the thermal conductivity. 

Eq. (3.17) may be considered as generalizing the 
method of the associated flow field of reference 4 to 
thermoelasticity by associating an entropy flow with a 
given stress configuration. 


(4) APPLICATION TO THERMAL STRESS ANALYSIS FOR 
STATIONARY TEMPERATURES 


Application of the reciprocity relations to thermal 
stress analysis may be made as follows: Consider a 
structure as shown in Fig. | and suppose we wanted to 
determine the stress at a plane cross section 7 due to 
the application of a temperature #p per unit area at 
point P of the boundary. As mentioned above, this 
temperature may be applied directly to the solid bound- 
ary or through a surface heat-transfer layer. The 
stress at / is represented by a normal force F,, a tan- 
gential force F;, and a moment SN. We cut the struc- 
ture along the plane 1/ and apply such forces very 
slowly—i.e., isothermally at the cut. The linear and 
angular displacements at the cut associated with the 
forces F,, F,, and are, respectively, u,, and a. 
They may be written 


Un = Cols CuFt + Cus 
= + CuFt + ( (4.1) 
@ = + + Cm 


The matrix of coefficients is, of course, symmetric, and 
relations (4.1) may be conveniently expressed by apply- 
ing Castigliano’s principle—i.e., 
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u, = OU /OF,, OU OF,, a = OU/OM (4.2) 


where is a quadratic form representing the 
isothermal strain energy. Up to now we have followed 
the well-known classical procedure. We now suppose 
that we have evaluated for each of the forces F,, F,, 
and IN, the entropy flow, and, in particular, the normal 
inward flow of entropy per unit area at point P. For 
each force component this entropy inflow at P may be 
written 


Sp! n) A PM F,, 
Sp' t) A DF; 


( (4.3) 
= A pa 


In evaluating the flow field S we must assume that it 
remains continuous across the gap at M. The calcula- 
tion is conveniently carried out by a variational pro- 
cedure as outlined in the previous section. We are 
now in a position to apply the reciprocity relations 
(3.6). This means that we may evaluate immediately 
the displacements produced at the gap by a temper- 
ature 6p applied at point P per unit area. We write 


“uy = A 6p, 
a= (4.4) 


un, = 6p, 


The influence coefficients in these relations are already 
determined by Eqs. (4.3) since from the reciprocity 
relations (3.6) we find 


= 


= Apy™ (4.5) 


Ayp” = Apyw™, 


Thermal stresses due to the temperature 6p correspond 
to forces F,,, F;, MM which, at the gap, produce displace- 
ments which are equal and opposite in sign to those 
Eqs. (4.4) produced by 6p. Hence the thermal stresses 
at ./ are given by the equations 

Apy6p = —(0U/OF,)) 


A —(0U/OF;,) > (4.6) 
= —(0U/OM). 


The reader will note the advantages of this procedure 
over the classical method. They are as follows: 

(a) The method does not require the knowledge of 
the temperature field in the body and by-passes com- 
pletely the necessity of calculating this temperature 
distribution. 

(b) The thermomechanical influence coefficients 
Apy, Apu, appearing in Eqs. (4.6) are 
determined for all points P of the boundary by a single 
calculation. The calculation does not have to be re- 
peated for every new distribution of boundary tem- 
perature. 

(ce) The thermomechanical influence coefficients are 
determined by a variational method which amounts 
to minimizing the dissipation function. This function 
includes the effect of any surface heat-transfer layer 
with a transfer coefficient which may be dependent 
on the location. 

(d) The method avoids the evaluation of compli- 
cated temperature fields which may arise due to local 
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effects near the points of application of the boundary 
temperatures, and provides a “smoothing out’’ of these 
effects. 

Needless to say, the method is applicable to very 
complex structures and is not at all restricted to the 
particular configuration of Fig. 1 which is used here to 
explain the method. As shown by this writer some 
years ago,’’* the connectivity of a body introduces 
essential differences in the thermoelastic properties. 
In problems of two-dimensional stress, for example, 
we have shown that an homogeneous isotropic body 
without holes—i.e., simply connected—exhibits no 
thermal stresses under arbitrary stationary temper- 
atures at the boundary. If there are holes, a cut such 
as exhibited in Fig. | will render the body simply 
connected and thermal stresses will disappear. They 
will then be due entirely to the forces required to close 
the gap at the cut. Application of these theorems to 
the photoelastic determination of thermal stresses was 
discussed in reference 9. 

The reciprocity relations are equivalent to a general- 
ization of Castigliano’s principle for thermomechanical 
phenomena, if we introduce the concept of associated 
field as developed in reference 4. The thermoelastic 
potential may be written 


V = (1/2) (u,’F, + + + Spbp) (4.7) 
where 


u,’ = OU OF, + Op 
= OU/OF, + Aur” Op 
a = 0U/OM + (4.8) 
Sp = F, + Apu + 
Apy@mM + App Op 


The quantity Sp or total entropy inflow at P is the sum 
of expression (4.3) and an additional term A pp§p which 
is the inflow due to @>p itself after elimination of the ig- 
norable coordinates. This elimination may be ac- 
complished by introducing associated flow fields as 
developed in reference 4. The quantities u,,’, u,’, and a’ 
are the total displacements at 7 obtained by adding 
Eqs. (4.2) and (4.4). Because of reciprocity and the 
properties of quadratic forms we may write 


u, = OV/OF,, 
a’ = OV/OM, 


u,’ = OV/OF,, 
Sp = (4.9) 


This is a generalized form of Castigliano’s principle. 
Putting 
(4.10) 


is equivalent to Eqs. (4.6) and determines the thermal 
stresses. Eqs. (4.10) also state that the thermoelastic 
potential is a minimum when considered as a function of 
a self-equilibrating stress field. 


(5) EXTENSION TO TRANSIENT STRESSES AND TO 
VISCOELASTICITY 


Until now we have only considered the application of 
reciprocity relations to the case of stationary temper- 
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atures. The next question is to extend the procedure 
to the analysis of thermal stresses in an elastic system 
for the case of transient temperatures. Since the ad- 
mittance matrix A,;* is symmetric, it is clear that a 
similar reciprocity property as Eq. (3.6) exists for 
transients. Considering again point P and M, and a 
temperature 6p(t) function of the time ¢ applied at P, 
we may write, for the displacement at 1/7, the opera- 
tional relation 


Uy (t) = A wp*Op(t) (5.1) 


In this expression A yp* is an operator of the type (3.2). 
This operator corresponds of course to an indicial 
cross-admittance function A yp(t). This admittance 
function is obtained as the displacement u, at .V cor- 
responding to the sudden application of a unit constant 


temperature at P. We put 


Op(t) = (5.2) 
where 1(¢) is the Heaviside step function 


 t<0 
I(t) = t> 0 (5.3) 


The indicial admittance is 
Awplt) = Axyp* 1(t) (5.4) 


Conversely, for the entropy inflow at P due to a force 
Fy,(t) at \/, we may also write the operational relation 

Sp(t) =x Ap y(t) (5.5) 
The corresponding indicial admittance A py(t) equal to 
the entropy inflow at P due to the sudden application 
of a unit force Fy at Vis 


Apu (t) = Apy* (5.6) 
The svmmetry of the admittance matrix (3.2) leads to 
Ayp* = (5.7) 


and, hence, to Axyplt) = Apy(t) (5.8) 


This is a reciprocity relation for the indicial thermo- 
mechanical admittance functions. 

Application of these results to thermal stress anal- 
ysis requires the evaluation of the indicial admittance 
Apy(t). This requires the evaluation of the transient 
entropy flow field due to the sudden application of a 
unit force at M. Calculation of this field may be 
achieved fairly simply by an approximate method. 
The general thermoelastic theory’: establishes a rela- 
tion between the temperature @ and the deformation, 
which is 


(0/dx,) = 
c(00/Ot) + (5.9) 


The deformation e,; is due to the suddenly applied unit 
force at W/. Since we are dealing with a thermoelastic 
system, if we wish to be exact the deformation e,; must 
be computed along with the coupled temperature field 
and will be a function of time. Actually, however, the 
coupling between the deformation and the temperature 
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is very small, and it is justified to assume that the de- 
formations e,; are the same as if they were isothermal. 
At this point, we make use of a scalar concept already 
introduced in reference 3, which we called a flow po- 
tential. In this case, we define it as 


¥ = (1/T,) oa (5.10) 
0 


Eq. (5.9) may then be written 
(0/0x;) [ki(OW/Ox;) | = c(Oy/Ot) + (5.11) 


Because of Eq. (3.16), we may also write 6,e;; = ajjoi; 
where o;,; is the stress due to Fy. 
The scalar y plays the role of a temperature field 


with distributed heat sinks Bye; = aijox per unit 
volume. The entropy displacement is then given by 
Si = (5.12) 


Actually it is not necessary to evaluate the scalar y, 
and the field S; may be computed directly by vari- 
ational methods as developed in considerable detail in 
references 3 and 4 in connection with purely thermal 
problems. The problem here is the same as the evalu- 
ation of the heat flow, represented by S when a sudden 
distribution of heat sinks 8,,e;; is applied throughout 
the body. This can be done by considering the prob- 
lem as one of thermal relaxation, by splitting S into 
two parts, 


S=S,+S, (5.13) 


where Sp is the steady-state flow corresponding to the 
sources, while S, is a transient field without sources 
determined by the initial condition S, = —S») and 
vanishing fort = ©. The steady-state field S») has 
already been considered above in connection with sta- 
tionary problems and, as is shown, may also be evalu- 
ated by variational methods. 

We now go back to the more specific configuration of 
Fig. 1. By the method just outlined we may evaluate 
the three admittance functions A yp” (t), 
Awyp'(t), which are associated with the normal tan- 
gential and angular displacement at the cut \/. A 
time varying temperature 6p(t) applied at P per unit 
area produces at M displacements given by 


u,(t) = Amp (t)6p(O) + | 
t 
(t T) (d0p/dr)dr 

0 


u,(t) A (t) @p(0) + 
f A yp'?(t — 1) (dOp/dr)dr 
0 
a(t) = (t) 6p(0) 
Aup™® (t — r) (dOp/dr)dr 
0 


| 
(5.14) 


| 

If, again, we neglect the temperature changes produced 
by the forces F,, F:, MN applied at M, the forces neces- 
sary to close the gap are the same as in Section (4) for 
isothermal deformation. Using the same isothermal 
strain energy U we write 


A yp (t)6p(0) + Ayp(t — 1) X | 
0 

(d0p/dr)dr = —(OU/OF,) 


| 


A (t)0p(O) + (t— X 
0 (5.15) 


| 


(d6p/dr)dr = —(OU/OF,) | 


| 

t | 

A + Amp'*(t — r) X | 
0 

(dOp/dr)dr = —(O0U/O0M) | 


These three equations yield the thermal stress at WM asa 
function of time. 

The reciprocity relations and the above methods are 
readily extended to what might be called linear thermo- 
viscoelasticity —i.e., the problem of thermal stresses in a 
viscoelastic material such that the operational moduli 
are independent of the temperature. In such a case the 
stress-strain relations are formally identical with Eqs. 
(2.4) where the coefficients c* , Bij, are replaced by the 
operators 8:;*. The operators and are 
related to the coefficient of thermal dilation by equa- 
tions similar to (3.14); i.e., 


Cay = (5.16) 


The requirement that the operators be temperature- 
independent strongly restricts the applicability of 
thermoviscoelasticity since most viscoelastic materials 
show deformation rates which are very sensitive to tem- 
perature changes. We shall therefore only briefly 
indicate how the method may be extended to this case. 

The validity of the reciprocity relations (5.8) for the 
thermomechanical admittance operator in viscoelastic 
media may be immediately derived from either the 
correspondence rule or the general thermodynamic 
principles as formulated in earlier publications. The 
admittance functions are again computed by solving 
Eq. (5.11) for Y or S. The sinks in that equation due 
to the application of the force Fy are 8;;*e;; = ajjoj,. 
We see that if o,; is statically determined they are con- 
stant. The admittance functions thus calculated are 
related to the stresses at MJ by equations formally 
identical with (5.15) except for the fact that the co- 
efficients on the right-hand side are operators. The 
reader familiar with operational or Laplace transform 
methods will have no difficulty in evaluating these 
stresses by standard procedures. 

Castigliano’s principle may be extended in oper- 
ational form to the transient thermoelastic case by 
writing an expression identical in form with (4.7) ex- 
cept that the coefficients are replaced by the corre- 
sponding operators. The same generalization holds for 
thermoviscoelasticity or any other more complex field 
of application of irreversible thermodynamics. It also 
corresponds to a generalization of methods of comple- 
mentary energy to thermoelasticity and thermovisco- 
elasticity. 


(6) EXAMPLE 


We shall consider the problem of stationary thermal 
stresses in a thin circular cylinder of radius a and thick- 
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ness h (see Fig. 2). The deformation is restrained along 
the axis, and the problem is one of two-dimensional 


strain. The two-dimensional stress-strain relations are 
Orr = + rE — 
= + Ae — (6.1) 
Ory = 
with the dilation e = e,, + e,, and Lamé constants ) 
and uw. The stress o.. in a direction parallel with the 
axis 1S 


0.2 = re — BO (6.2) 


Following the procedure outlined above, we cut the 
cylinder at 7 (Fig. 2). Forces and moments applied 
at 7 produce a bending moment Mp at point P. 


Mp = M+ F,a(l — cos ¢g) + Fasin ¢ (6.3) 


The forces are expressed per unit length along the axis. 
The strain energy due to these forces is 


U = (1/2)[0 — »°)/E| | 


0 (6.4) 
U = (1/2)R[290? + 3a°F,? + + 


with 


| 


= [(1 — v*)/E| (12/h') za 
E = Young’s modulus 
v = Poisson’s ratio 


Castigliano’s principle expressed by Eqs. (4.2) yields 
the linear and angular displacements at M. 


u, = 3a°RF, + 2aRM| 
= Ra*F, (6.5) 
a = 2aRF, + 2RM ( 


We must now calculate the displacements at WM due to 
temperatures at the cylindrical boundary by applying 
the thermomechanical reciprocity relations. The forces 
at M produce a dilation e and a corresponding entropy 
flow S which we assume to be perpendicular to the 


circumference. The region around point P is repre- 
sented in Fig. 3. The flow S across the thickness must 
satisfy 


0S/dy = —Be (6.6) 


where the dilation ¢ is linearly distributed and given byt 


= L 7) 

D = [(1 + ») (1 — 2v)/E] (12/h3) 
Integrating Eq. (6.6), 

S = (1/2)eDy? + C (6.8) 


where C is a constant of integration to be determined 
by a principle of minimum dissipation—i.e., by mini- 
mizing expression (3.23). If there is no heat-transfer 
surface layer (K = ~) this results in the equation 


+(h/2) 
(0/dC) Stdy = 0 (6.9) 


— (h/2) 


+ We shall neglect here that part of the deformation which is 
due to other forces than the bending moment Wp. 


y 


Fic. 2. Stresses in a thin circular cylinder. 


Fic. 3. Section across the thickness at P. 


from which we determine C. We find 
S = (1/2)8D[y? — (6?/12)} 


The inward component Sp of S at the inside boundary 
(y = —h/2) is 


(6.10) 


Sp = (1/12)BDh* (6.11) 
Sp = + v) (1 — 2v)/EA\ Mp 
It is easily shown that 
y = (@/E) (1 + v) Cl — 2p) (6.12) 


represents the coefficient of thermal dilation for two- 
dimensional strain; hence, we may write 


Sp = y(Mp/h) (6.13) 


For each force component at / we derive from Eq. (6.3) 


Sp” = (ya/h) (1 — cos Fa) 
Sp (ya/h) F, sin ¢ 
Sp = (y/h)m 


These equations correspond to relations (4.3). The 
thermomechanical influence coefficients are 


A yp” (ya/h) (1 — cos 
Amp” = (ya/h) sin ¢ 
= y/h 


We shall apply relations (4.6) to the case of an arbi- 
trary distribution of temperatures @p along the inside 
circumference. On the left-hand side we must there- 
fore integrate the temperatures over all points P. 
The right-hand side is given by relations (6.5). We 
find 


(6.14) 


(6.15) 


Aup™0pade = 3a?RF, + 2aRM 
AypOpade = Ra®F, (6.16) 
= 2aRF, + 2RM 

If the temperature is uniform, i.e., 0p = 6%, we find 
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= () and 
M = [E/(1 — v?)|p(h*/12) (6.17) 


The thermal stress is a pure moment. For 6p = 
4 sin yg, we find F,, = M = Oand 


F, = —(76/ah) (E/(1 — v?)] (A8/12) (6.18) 


For 6p proportional to sin ng or cos ng where (” = 
integer) 7 > 1 the thermal stress vanishes. These re- 
sults are in agreement with our general theorem of 
reference 8 on two-dimensional thermal stresses. 

If, instead of applying the temperature 6p directly 
to the solid boundary, we apply it through a layer of 
heat-transfer coefficient K we must minimize the ex- 
pression 

+(h/2) 


D’ = (1/2k) S? dy + (1/K)Sp? (6.19) 


—(h/2 
with 


k = heat conductivity of the solid 
S = (1/2)6Dy? + C (6.20) 
Sp = (1/8)BDh? + C ( 


The factor 1/2 does not appear in the second term of 
Eq. (6.19) because it is the sum of two terms corre- 
sponding to the inner and outer boundaries. _Minimiz- 
ing D’ determines C, and we find 


Sp = BDh?/12[1 + (2k/Kh)] (6.21) 


Comparing with (6.11), we see that all calculations may 
be repeated, provided we multiply all thermomechanical 
coefficients A yp", ete., by 1/[1 + (2k/Kh)]. 


JOURNAL OF THE AERO/SPACE SCIENCES—JULY, 1959 


The example presented here has been chosen for its 
simplicity and as an illustration of the methods. It is 
not intended to show the particular advantages attached 
to the procedure ‘This should come out in the treat- 
ment of more complex problems. 
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A Theory of the Response of Airplanes to 
Random Atmospheric Turbulence 


B. ETKIN* 


Unwersity of Toronto 


SUMMARY 


A theory is given for the aerodynamic forces and moments 
which act on an airplane as a result of flight through turbulent 
air, and for the associated dynamic response of the airplane. It 
starts with Ribner’s spectral theory of the phenomenon, and ex- 
tends it to a specific formulation of the forces and moments. This 
is achieved by expanding the downwash field of a single spectral 
component in a power series in the neighborhood of the airplane. 
Spanwise variations of the gust intensity are included, and the 
wing can be treated by any applicable lifting line or surface 
theory. The final expressions give the forces and moments in 
terms of a set of aerodynamic derivatives like those used in sta- 
bility and flutter analysis. 

Consideration of the equations of motion, and the airplane 
transfer functions, shows that the final response quantities of 
interest can be obtained in the form of one-dimensional power 


spectra. Tables and graphs of the relevant input spectra are 
given. 
SYMBOLS 
A = rolling moment of inertia 
R = aspect ratio, b/c 
B = pitching moment of inertia 
b = wing span 
e = yawing moment of inertia 
( = wing reference chord 
e = equivalent elastic motion due to gust, 
= (c?/8u9) = —(1/2u9)ki2w,(0, 0) 
es = equivalent elastic motion due to gust, 
= (bc/4uy) (07w,/Ox = —(1/uo)Rikow, 0, 
e = equivalent elastic motion due to gust, 
= (b?/8u9) (0?w,/Oy?)o = —(1/2uo )Ro2w,(0, 0) 
E = product of inertia /,, 
fi, fo = slope correction factors [Eq. (2.8)] 
i4 = 84/pSb' 
iz = 8B/pSc* 
= 8C/pSb* 
ip = S8E/pSb* 
le = generalized inertia in elastic mode 
ky = nondimensional wave number 2\c/2 
ko = nondimensional wave number 22.)/2 
k. = nondimensional oscillation frequency of elastic 
mode 
l; = length of wing (Fig. 2) 
Lb = scale of atmospheric turbulence 
m = airplane mass 
Pp = roll rate, rad./sec. 
Po = equivalent roll rate due to spanwise gust gradient, 


= —fo( OW, )o 
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Do = pgb/2uy = 0, 
q = pitch rate, rad./sec 
Go = equivalent pitch rate due to streamwise gust 


gradient = /\(Ow,/Ox)o 
qo = gGyl/2uy = iki fiw, (0, 


r = yaw rate, rad. /sec. 

= wing area 

t = time 

i = 2uct/c in longitudinal equations; 2%of/b in lateral 
equations 

Uo = airplane forward speed 

Wg = gust velocity (positive downwards) 

w (0,0) = value of w, at airplane center of gravity 

M1, Xe = space-fixed coordinates 

= moving coordinates 

i = 2x/c 

AY = 2y/b 

Oe = equivalent angle of attack due to gust, 
—w,(0, 

= sideslip angle 

o = root mean square value of w,, ft. per sec. (gust 
intensity ) 

€ = generalized coordinate of elastic mode 

ny = wavelength of spectral component (Fig. 1) 

Mt, Az = component wave lengths (Fig. 1) 

Q = wave number vector 

Q, Qe = components of Q 

Q', 22’ = cut-off wave numbers (see Fig. 11) 

¢ = angle of bank 

$(2), 22) = two-dimensional power spectrum of 7w, 

(21) 

2/2) = derived one-dimensional input spectra {see Eqs. 

(8.3) and (8.4)] 

) 

m = 2m/pSc in longitudinal equations, 2m/pSb in 
lateral equations 

p = air density 

¢ 3 = indicates Laplace transform with respect to / 


(1) INTRODUCTION 


Tt PROBLEM of calculating the response of an air- 
plane to atmospheric turbulence contains several 
distinct elements. These are (1) the statistical de- 
scription of the turbulent field (the input), (2) the 
calculation of the aerodynamic forces associated with 
the turbulent field (the ‘‘gust forces’’), (3) the calcula- 
tion of the transfer functions which relate the airplane 
motion (or other quantity of interest, such as stress) 
to the gust forces, and (4) the combination of the input 
with the transfer functions to obtain the output (mo- 
tion, stress, etc.). 

The details of the steps (2) to (4) are strongly in- 
fluenced by the method adopted for describing the 
turbulent field in step (1). Two basically different 
approaches have been taken to this description. The 
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Fic. 1. The elementary spectral component. 


first was introduced by Liepmann' * and was subse- 
quently developed by Diederich® and others. It 
utilizes as its basic element the two-point correlation 
function of the velocity component of interest. The 
second approach was introduced by Ribner.‘ It repre- 
sents the turbulence by a superposition (Fourier in- 
tegral) of sinusoidal waves of shearing motion. 

This paper presents an analysis which utilizes Rib- 
ner’s formulation for the turbulence. It deals only 
with the vertical component of the gust velocity (wy, 
positive downwards). 

The method is based upon the approximation of a 
single spectral component of the turbulent downwash 
by a modified two-dimensional Taylor series in the 
neighborhood of the airplane. It is found that only 
terms up to the second degree need be retained in order 
to obtain a good fit for wavelengths down to twice the 
relevant airplane dimension (span or length of wings). 
Truncation of the spectrum of the turbulence at this 
wavelength, and the discarding of all the components 
of shorter wavelength, reduces the root mean square 
downwash by a small amount (typically about 2 to 3 
percent). This is negligible in view of the uncertainty 
in the data on atmospheric turbulence. The omission 
of the short-wavelength components of the turbulence 
also results in cutting off the spectrum of the input to 
the airplane above a certain frequency. This cut-off 
frequency is found to lie above the frequency of the 
lowest structural modes of high aspect ratio airplanes, 
for which these modes are important. Thus, the 
response of the airplane in its structural modes is in- 
cluded within the scope of the analysis for those cases 
in which this response is important. 

The aerodynamic coefficients appear in the equations 
in the form of transfer functions (or frequency-response 
vectors), of precisely the same kind as are commonly 
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used in stability and flutter analysis. Thus the 
method presented is a natural extension of existing and 
well-developed techniques. 


(2) THE REPRESENTATION OF THE TURBULENCE 


Summary of Previous Results Required for the Analysis 


As noted above, the representation of the turbulent 
field is that due to Ribner.* In it, the downwash pat- 
tern in the x, x2 plane is obtained by the superposition 
of infinitely many elementary inclined waves of shear- 
ing motion (see Fig. 1). These elementary waves are 
of all orientations, and all wavelengths. Each is given 
by the expression 


i(Qix1+Qex2) 


dw, (x1, Xe) = ae —o < Qi, < @ (2.1) 


in which a is the (infinitesimal) amplitude of the wave, 
and 


Q. = 2r de (2.2) 


are the wave number components and )j, A» are the 
wavelengths on the two axes (see Fig. 1). It should 
be noted that the usual assumption of a ‘‘frozen’’ field 
is adopted here. That is, the local downwash varies 
with position, but not with time. 

The contribution of the spectral component to the 
mean square downwash (averaged over the x, x2 plane) 
is given by 


dw,? = $(Q:, %)d%dQ, (2.3) 


where $(Q), Q2) is the two-dimensional power spectrum 
of the turbulence. Thus, the mean square downwash 
from all the spectral components is 


v2 =o = f $(Q), 22)dQdQ, (2.4) 


The available information on atmospheric turbulence 
indicates that it may be approximated as homogeneous 
and isotropic over regions of limited extent.’ The 
spectrum has moreover been found to be described 
adequately by the same function that applies to wind- 
tunnel turbulence. The two-dimensional spectrum 
function is 
22) = (3/4) (o?/r)L? (LQ)? X 

[1 + (LQ)? (2.5) 
where Q? = + 


and ZL is the scale of the turbulence. Integration of 
the above expression with respect to 2 provides the 
one-dimensional spectrum function—-i.e., 


$1(21) f 


(o?L/2m) { [1 + 3(L)2] + 
[1 + 


(2.6) 


Press and Meadows? suggest a value of L = 1,000 ft. 
as representative of atmospheric turbulence, and also 
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give some data on the probability distributions of o 
for various operating and weather conditions. 


The Series Approximation to the Spectral Component 


The basic problem is the calculation of the aero- 
dynamic forces and moments acting on an airplane 
when flying through an inclined sinusoidal wave of 
downwash having unit amplitude—ie, w, = 
+ To carry out this calculation in an 
exact form for airplanes of arbitrary plan form and for 
all Mach Numbers seems to be a formidable task. 
In particular cases, numerical integration of the lifting 
surface integral equation can provide the required 
aerodynamic transfer functions. In this paper an 
approximation to the forces and moments is obtained 
which is based on a series representation of the down- 
wash in the neighborhood of the airplane. The series 
used is 


w(x, V) = (Wy)o + fi(Ow,/Ox)ox + 
fo(Ow,/Oy)oy + (1/2) (0? w,/Ox*) ox? + 
(0°w,/Oxdy)oxy + (1/2) (A*w,/Oy*)oy? (2.7) 


Here the subscript zero indicates that the derivatives 
are evaluated at the origin of the xy system-—i.e., at 
the airplane center of gravity. The factors f; and fe 
are introduced to improve the fit obtained with Eq. 
(2.7) [see Section (4)}. They are 


fy = 1 — 0.084 (Q)/,/2)? 


where /; and /, are the lengths shown in Fig. 2. 
The coefficients of the power series are computed as 
follows :* 


(2.9) 


#(Qix + Qey 
hence, w, =e 2.10) 


The coefficients of Eq. (2.7) follow directly as 


(Ow,/Ox)o = iQ,w,(0, 0) 
(Qw,/Oy)o = i2.w,(0, 0) (2.11) 
(0?w,/Ox*)) = —Q,? w,(0, 0) | 
(0?w,/Oxdy)o = —MQw,(0, 0) 
(0?w,/Oy")o = —2* w,(0, 0) 

and hence Eq. (2.7) becomes 

w(x, y, 2) = [1 + ifiQix + ify — 

(1/2)Q;2x? — QQeaxry — (1/2)Q*y?} (2.12) 


Thus, the downwash field is seen to consist of the super- 
position of several relatively simple spatial distribu- 
tions, each of which varies periodically with circular 
frequency (4%. Eq. (2.7) may also be rewritten in 
terms of the nondimensional quantities defined in the 
list of symbols as 


* The statistical properties of the input to the airplane are 
obtained by assuming that the airplane has a rectilinear motion 
at constant speed through the gust field. 


4, => 


Fic. 2 


Wy = Uo(—ay + — Pov + Erk? + + 
(2.13) 


Here ay, Gy . . . €3 are periodic functions of time only. 
(3) INTERPRETATION OF THE TERMS IN THE SERIES 
APPROXIMATION 


It can be shown that when a flat-plate wing lying 
in the xy plane flies through a weak turbulent velocity 
field (frozen field) in which the turbulent velocity nor- 
mal to the wing is w = f(x, y, ¢), then the pressure dis- 
tribution on the wing is identical with that which 
exists when the same wing flies through air at rest, but 
has a surface slope and velocity such that the boundary 
condition is w = —f(x, y, 2). 

This boundary condition, together with the appro- 
priate differential equation, completely defines the 
mathematical problem of calculating the wing load dis- 
tribution. The total load distribution corresponding to 
the boundary condition represented by Eq. (2.13) is 
seen to consist of a sum of loadings, one for each term of 
the equation. The boundary conditions corresponding 
to the various terms are individually of quite simple 
form. 

Each term of the expression for the gust velocity 
[Eqs. (2.12), (2.13) ] can be interpreted, from the stand- 
point of pressure distribution, as describing a certain 
periodic rigid or elastic motion of the aerodynamic 
surfaces. These modes of motion are discussed below, 
simply as a convenience in attaching a physical sig- 
nificance to the terms of the equation. 


bed] 

a, = —w,(0, 0)/m. The term w,(0, 0) of Eq. 
(2.12) represents a constant contribution to w,(x, y) 
and, hence, a uniform downwash over the xy plane. 
The equivalent motion is therefore an upward transla- 
tion of the airplane, and the loading is that due to the 
angle of attack a,. 


A 


Go 

The term = represents a downwash propor- 
tional tox. The surface motion producing the appro- 
priate boundary condition is a pitching velocity g, (see 
Fig. 3). Hence, the loading associated with this term 
is that due to an angular velocity in pitch of nondimen- 
sional magnitude g,. The 7 multiplying this term in 
Eq. (2.12) simply indicates that the equivalent pitch- 
ing is 90° out of phase with the angle of attack. 
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Bo in pitch of a magnitude proportional to j. This load- 


The term which is linear in y is analogous to the one 
linear in x. The equivalent surface motion is a rate 
of roll, of nondimensional magnitude /,. 


ey 

The term. wet” describes a parabolic downwash. 
The equivalent surface motion is a streamwise elastic 
flapping, as illustrated in Fig. 4. The periodic surface 
motion, in order to produce the required boundary 
condition, must satisfy the equation 


(O2/Ot) — Uo(OZ./Ox) = — (3.1) 


where 2,,.(x, ¢) is the instantaneous surface position, and 
é; is the periodic function 


= 
Let u(x, t) = 
then Eq. (3.1) becomes 
(dZ,/dx) — = (3.2) 
This equation has the solution, for Z,,(0) = 0, 


Z (xX) = (1/7Q\u) x 
— 1) — ix + (1/2) (Qix)?] (3.3) 


For small values of 2, this expression can be approxi- 
mated by expanding the exponential term to the third 
degree with the result 


Zeolx) = —(2x3/6u0) (3.4) 


The equivalent motion in question may thus be re- 
garded as a periodic cambering of the wing according 
to the relation (3.3), or at low frequencies, the approxi- 
mate equation (3.4). 


€2 

The term wmeé.xy describes a downwash varying 
linearly with both x and y. The equivalent surface 
motion is a periodic cambering at each section, as illus- 
trated in Fig. 5. The amplitude of the camber is pro- 
portional to y, and antisymmetric with respect to the 
centerline. The camber line equation at any section y 
is obtained from the boundary condition [e.f. Eq. 


(3.1) ] 
where Co = 


Solution of Eq. (8.5) gives the surface displacement 
function 


Zw(x, Y) = — 1 — iQyx) 


For small values of 9; this is approximately 


Zw = yx? 


When strip theory is being used to calculate the 
aerodynamic coefficients, the contribution of the e. mode 
becomes very simple, for then the loading at each sec- 
tion would be simply that due to an angular velocity 


ing would be known from two-dimensional theory. 


e3 

The term mesy? of Eq. (2.13) describes a downwash 
varying parabolically in the spanwise direction, and 
constant in the stream direction. The equivalent 
surface motion is a periodic spanwise flapping, as shown 
in Fig. 6. The shape of the surface is found as in the 
preceding cases, and is given by 


Lw = Qe?y?/ 212 (3.6) 


(4) AccuRACY OF THE APPROXIMATE DOWNWASH 
EQUATION AND LIMITATION ON THE WAVE LENGTH 


Any vertical cut through the exact sinusoidal down- 
wash field gives a sine wave for w, The same cut 
through the approximate field given by Eq. (2.7) gives 
a parabola. The approximation can therefore be re- 
garded as that of fitting a parabola to a sine wave. 
We examine the fit obtained along the x axis as being 
representative of any diagonal axis. The following real 
expressions then give the exact and approximate 
downwash functions, where the subscript zero denotes 
x = 0: 


Exact: W, = sin Qxy = sin Qy(uot + x) 
(Wy)o + fi(Ow,/Ox)ox + (4.1) 
(1/2) (07w,/2x*)x? 


Approx: W, 


or 


Exact: wy, = sin Qyot cos Qx + 
cos sin 
Approx: w, = sin — (1/2)Q)2x?] + 


MYX Cos 


(4.2) 


A number of tests are applied to the approximate 
function to fix the range of validity; these are described 
below. The conclusion arrived at is that the approxi- 
mate equation is satisfactory down to wavelengths 
equal to twice the wing length. The approximation 
for this extreme condition is shown in Fig. 7. 


Root Mean Square Value of the Approximate Downwash 


Let the interval over which the approximate function 
is applied be 


—(l,/2) < x 
Then the mean square value of w,, averaged over the 
above range of x, and over time, is 


1/2 
w,? = (1/2mh) ao dx w,?(x, 6) (4.3) 
0 


where 0 = Qyuot 


The result of the integration, when the approximate 
expression for w, is used, is 


w,? = (1/2) — (1/6) (Q,/2)? (A — fi?) + 
(1/40) 4 (4.4) 
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Fic. 9. Accuracy of moment of approximate downwash. 


This expression will give the exact value 1/2 if 
fy? = 1 — (4.5) 


The actual value used for f; (Eq. 2.8), is chosen to give 
a good approximation to this over the range of interest. 
The ratio of the approximate and exact root mean 
square value of w, is shown in Fig. 8. This graph shows 
that the average value of the gust intensity is given 
with negligible error by Eq. (2.13) for gust wavelengths 
down to twice 


Error in the Moment of wg 


The result above shows that, on the average, the 
approximate downwash is close to the exact downwash. 
It does not, however, constitute a test for the accuracy 
of the spatial distribution. In particular, the aero- 
dynamic rolling and pitching moments due to the gust 
are expected to be closely related to the corresponding 
moments of the w, distribution. The exact and ap- 
proximate values of this moment are therefore com- 
pared. The moment of the w, distribution about the 


origin is 
11/2 
M = f wx dx 
—1,/2 


When the two expressions for w, are used [Eq. (4.2) | 
we obtain after integration 


M approx. Ai 
3 sin(Q,/;/2) = (Qi /2) cos /2 

independent of ¢. This ratio is shown on Fig. 9, from 

which the maximum error in the moment is seen to be 


less than 2'/, per cent. Comparably small errors have 
been found in the ‘‘one-side’’ moment 


/ 2 
wx dx 
0 


to which all terms in w, contribute. (Note that only 
the odd terms in w, contribute to .) 


(4.6) 


The Root Mean Square Error 


A much more severe test of the approximation is 
obtained by considering the root mean square error. 
For this purpose, we consider the difference between 
the exact and approximate downwash 


e(x, t) = (Wp (Wy) exact 
sin — (1/2)Q)2x? — cos Qx] + 


cos Quot(fiQix — sin Qx) (4.7) 


The mean square error is obtained by averaging over 
the same space-time domain as used in Section (4); ice., 


11/2 
= (1/2mh) f dé f dx x) (4.8) 
0 


—h/2 
The result of carrying out the integration is 


= 1— [(1 — + (1/40)y4 — 
((2 + fi)/¥] sin + (1/2) sin + 
(1+ fi) cosy (4.9) 


where y = (Q/,/2). The root mean square error has 
been calculated from this equation and is plotted in 
Fig. 10 as a percentage of the root mean square value of 
w,. The abscissa is the ratio of fitting interval to wave- 
length. It is seen that the root mean square error 
amounts to about 12 per cent of o at the shortest wave- 
length used (/;/A; = 1/2). 

This error appears to be rather large. However, its 
effect upon the final results may be expected to be 
much less than 12 per cent. The error is this large 
only at the highest wave number, at which the gust 
intensity is lower by a factor of 10° (even for large 
airplanes) than at low frequencies where the error is 
very small. 


Reduction in Root Mean Square Downwash Because of 

Truncation of the Spectrum 

Figs. 8 to 10 indicate that the error involved in 
approximating the exact sinusoidal downwash by the 
series expansion used is acceptable for the spectral 
components of wavelength \; 2 2/;. The same cri- 
terion is applied in the spanwise direction—-i.e., \» 2 
2/2 = 2b. Thus in the plane of the wave number vector 
(see Fig. 11), all spectral components with wave num- 
ber vector lying outside the region S are excluded. 
This entails an error in the root mean square value of 
the downwash. The exact value, o, is given by the 
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Fic. 10. Root mean square error in downwash. 
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integral of Eq. (2.4). The approximate value of the 
mean square downwash is 


2,’ 
f 6(Q), 
—Q2’ 


The error involved is calculated approximately from 


f (21, Q2)dQ dQ, 


where .S’ is the region outside the circle C (Fig. 11), the 
area within which is the same as that of S. From Eq. 
(2.5), the expression (4.11) may be written as 


3 0 x (LQ)? 
L? 2 OdQ 
J. [1 + 


The range of interest is (LQ’) 2 5 for which the de- 
nominator of the integrand is to a close approximation 
(LQ)*. The integral has been evaluated with this ap- 
proximation, with the results shown in Fig. 12. A 
typical value of the error introduced is seen to be about 
2'/5 per cent. 


(4.10) 


(4.11) 


(4.12) 


The Cutoff Frequency 

The limitation on wavelength appears also as a cut- 
off in the input frequency [see Eq. (2.10)]. The limit- 
ing value is 


= tuo = Uo/ 21; (4.13) 


Fig. 13 shows the cutoff frequency as a function of 
size (/;), and speed. Also shown are three points for 
the wing-bending modes of representative high aspect 
ratio airplanes. It is for this class of airplanes that 
wing-bending flexibility must be taken into account in 
calculating the stresses due to gust loading. The cut- 
off frequency is seen to be higher than the lowest wing- 
bending frequency for all significant speeds. Thus the 
corresponding elastic modes can be included within 
the scope of the method. 


(5) THe Input Forces AND MOMENTS 


Associated with each term of Eq. (2.13) is a periodic 
pressure distribution acting over the surface of the air- 
plane. These integrate to give certain forces and 
moments. The relationships involved are conveni- 
ently depicted in the block diagram below (note 
that the X force is not included, since only five 
of the rigid-body degrees of freedom are retained). 
In this diagram C; (¢) is an appropriate generalized 
force coefficient for the wing bending mode (or other 
elastic mode). 
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The force and moment outputs are related to the in- 
put variables, a, .. . ¢;, by a set of transfer functions or 
complex stability and flutter derivatives*—i.e., 


C,(s) Gyp Pal 


Ci(s) = Gip Gres (5.1) 
C,(s) Gnp 
C.(s) & & Taw 
Cais Gus Ge, Gon (5.2) 
és (s) 


Eqs. (5.1) give the forces for the lateral equations of 
motion, and Eqs. (5.2) those for the longitudinal equa- 
tions. The omission of a number of terms; e.g., Gng, 
G.,,; is on the basis of symmetry. Only the symmetric 
(asymmetric) members of the input produce symmetric 
(asymmetric) forces and moments. The bar signifies 
the Laplace transform of the indicated variable-—e.g., 


0 


where 7 is the nondimensional time. 

A number of the transfer functions in the above 
equations are just those which ordinarily appear in 
the rigid-body equations of motion (see reference 6, 
Eqs. 11). Methods for determining these are dis- 
cussed in reference 6. The remaining transfer func- 
tions can be determined in the same way from frequency 
response measurements or theoretical calculations. 
For example, to obtain G,,,, the information required 
is the amplitude and phase of the Z force (the negative 
of the lift) of the wing flapping periodically in the mode 
e; (see Fig. 6). The calculation of such quantities is a 
common task in flutter analysis; the same techniques 
are applicable to the present problem. 


Contribution of a Tail 


Significant contributions, particularly to the pitching 
moment, may be expected to come from the tail, when 
there is one. The downwash at the tail consists of two 
parts, that due to the gust (w,);,,i; and that from the 
wing (W@w»)tai- We introduce an auxiliary reference 
system for the tail, as in Fig. 14. The gust downwash 
distribution at the tail for the unit spectral component 
may then be written 

Comparison of Eq. (5.4) with Eq. (2.10) shows that the 
direct effect of the gust on the tail may be treated in the 
same way as the effect on the wing, except that a phase 
lag of Q,/, appears. Since the tail is ordinarily con- 
siderably smaller than the wing, it is likely that in 
many cases only the (a@,)tai; distribution will need to 
be included. 

The exact determination of the wing downwash field 
at the tail presents a formidable problem. It is again 
reasonable, because of the smaller size of the tail, to 
assume a uniform downwash over it equal to the value 
at the mean aerodynamic center. This downwash will 
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be periodic, with the same frequency as that of the gust, 
and will have contributions from each of the symmetric 
wing distributions ¢:, and e3. For the low fre- 
quency components only the a, contribution is impor- 
tant. However, as the frequency increases, the others 
become relatively more important. Information on 
the periodic downwash in the wake of a wing seems to 
be scarce. This would seem to be a subject for future 
theoretical and experimental research. 


(6) Tue LoncirupINAL EQUATIONS AND TRANSFER 
FUNCTIONS 


The Laplace transforms of the longitudinal equations 
of motion, including one symmetric elastic mode (gen- 
eralized coordinate ¢) are written out in matrix form 
below, with the customary linearizing assumptions, and 
with the speed assumed constant at mw. On the right- 
hand side, the expressions (5.2) have been converted 
into a form in which w,(0, 0) appears as the input, by 
means of the definitions of a,, ete., and Eq. (2.11); 


a 
| = {Bt w,(0, 0) + w,(0, 0) (6.1) 
€ 
where = 
(Qus — Gig) —(2u + Giz) —G,, 
— (ips Ging) 
+ thiftGeg — | 
(BY =| [—Gna + ikifiGng — (1/2)R1?Gine,| | (6.2) 
[—Gsq + | 
1 
ic} 9 
Ge, 


The functions required for calculating the response 
to atmospheric turbulence are the squares of the moduli 
of the overall transfer functions, evaluated for periodic 
inputs i.e., with s = 7k;. These are 


= 9 9 


@,(0,0)), @,(0,0) , @,(0, 0) 


We consider now the manner in which these depend on 
the two nondimensional wave number components ; 
and k,. On the left-hand side of Eqs. (6.1), k, does not 
appear, neither explicitly nor implicitly. However k; 
appears explicitly (s = 7k,), and in general the aero- 
dynamic transfer functions are frequency dependent— 
€.g., Ging = Gma(tki). On the right-hand side of Eqs. 
(6.1), and both appear explicitly, and, moreover, 


f; and the aerodynamic transfer functions are also func- 


tions of k;. It follows then that the transfer functions 


are given by 


a/@,(0, 0) 
(0, 0) | = + c(R,)} (6.3) 
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where the notation indicates that the elements of the 
column matices on the right-hand side are functions of 
f;, but not of k2. These matrices are given by 


= (1/uo) [A] 


(6.4) 


_ Finally, the squares of the moduli of the transfer func- 


tions, which are the quantities needed for the response, 
can be calculated from Eqs. (6.3) and (6.4). That 
for the angle of attack, for example, has the form 


AIRPLANES TO 


ATMOSPHERIC TURBULENCE po 
a 
@,(0, 0) a’(k;) + 4 (ky) (6.5) 


The functions a’(k;), 6’(ki), c’(ki) can be calculated 
from Eggs. (6.3) and (6.4) for any given case and plotted 
for each of the response variables a, q, and e. 


(7) THE LATERAL EQUATIONS AND TRANSFER 
FUNCTIONS 


By proceeding as in the previous section, we obtain 
the lateral equations of motion in the form 


B 
uo [A’]| & | = ke {B’} w,(0, 0) + ko® {C’} w,(0, 0) (7.1) 
? 
(2us — — (Gos t+ Cr) (28 — Gr) ]) 
where [A’] = — Gig (i4s? — Gis) — (igs + Gy) 
— (igs? + Gyps) (tics — Grr) | 
IGyp — RiGye, | 
IGnp | 
ic} = —0. 084i [ Gy» | | 
| 
} 


(The cubic term in k, occurs because /, contains k,*.) 

Since the nondimensional time scale used in the lateral 
equations is based on the span, not the chord, then the 
input quantity w, is here given by 


__ 


w,(0, = ¢ 


The required frequency response functions B/@,(0, 0), 
etce., are calculated from Eq. (7.1) by replacing s with 


iki R. The transfer relation paralleling Eq. (6.5) 
for a typical variable-—e.g., 8-—is 
= (ki) + + 

(1) (7.3) 


Again, the functions a”(ki), 6”(ki) and can be 
calculated from Eqs. (7.1) and (7.2). 
(8S) THE MEAN SQUARE RESPONSE 


A typical mean square response in one of the longi- 
tudinal degrees of freedom -e.g., a—is of the form 


= 

$2(Q)) 

= [1 + — 22(Q)}} 
where 


(3/2) + (£9)? + (1/3)g%(2)) | 
(1/2) (o? {(L%)2/[1 + + (3/4) (02 + g(Q) 1/1 — g(Q)]} 


115 + log. { [1 + — g(Q)]} + (2/3)g%(Q) + 2[2 + (LM)? 
g(Q:) = (LQ")/[1 + (LQ)? + 


12)’ Qo’ 
a’ = f la @,(0, 0) (8.1) 
— i)’ —Q2’ 


where ¢(), Q2) is the power spectrum of the turbulence 
given in Section (2). 

In view of Eq. (6.5) and the definitions of k; and ky, 
Eq. (8.1) can be written 


(8.2) 
Qo’ 
(O91, Qe)dQ. 
— Qe’ 


—Q2’ 


where 


Qe! 
— 


Using Eq. (2.5) for Q), the integrals for . . . 
give 


2g(M%) — (2/3)g*(Q) (8.3) 


These one-dimensional spectrum functions have been calculated for three ratios of span to turbulence scale and are 


plotted on Figs. 15 to 17 and tabulated in Tables | to 3. 
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The mean square response in a lateral degree of freedom is calculated in a similar way. 


However, Eq. (7.3) con- 


tains a term in k,®, so that we require for the lateral response calculation a fourth one-dimensional spectrum, 


Q2’ 
= f 
— 2h’ 


Integration of this expression yields the result 


0.6 Q:, Qo) d Qe 


= (LO)? [1 + (} — — (5/2) log. [1 + — g(@)]} + 
+ (2/3)g*(M)) + (Bo? [1 + g(Q1)/2[1 — 12} — (13/4) — ]} + 
(35/8) log. {1 + g(%)]/[1 — g(Q)]} — 6g(M%) — (2/3)g%(M)) (8.4) 


This function is plotted on Fig. 18 and tabulated in 
Table 4. 


(9) A Test OF THE APPROXIMATION IN 
Two-DIMENSIONAL INCOMPRESSIBLE FLOW 


The approximate formulation for the aerodynamic 
forces is tested by calculating C, and C,, for a flat-plate 
two-dimensional airfoil in incompressible flow (the 
moment axis is at the mid-chord). The exact solution 
for this case has been given by Sears.?’ The approxi- 
mate solution in accordance with the present method is 
{see Eq. (5.2) ] 
= + + Greer \ 


; (9.1 


Employing the definitions of a,, gy, ¢:, these may be re- 


written as 


(0, 0) = —Geg + ikifiGeg — 


uo Cn 0) = + ftGing — 
TABLE 1 
Spectrum Function ¢;(Q) for L = 1,000, 7 = 10 
LQ, Lo’ = 10 = 50 = 100 
0.01 1,568 1,591 1,592 
0.02 1,569 1,591 1,592 
0.05 1,572 1,595 1,595 
0.1 1,583 1,606 1,607 
0.2 1,624 1,647 1,648 
0.5 1,759 1,782 1,782 
1 1,568 1,591 1,591 
2 805 827 827 
5 159 178 179 
10 33.1 46.0 46.7 
20 11.08 11.66 
50 0.374 1.349 1.707 
100 0.0475 0.213 0.388 
TABLE 2 
Spectrum Function @2(%) for L = 1,000, 6? = 10 
LY LQ.’ = 10 = 50 = 100 
0.01 8.00 XK 107° 1.562 * 107? 1.893 X 107? 
0.02 8.00 1.562 1.893 
0.05 7.99 1.562 1.893 
0.1 7.99 1.562 1.892 
02 7.97 1.559 1.890 
0.5 7.79 1.541 1.872 
1 1.477 1.807 
2 5.56 1.306 1.636 
5 2.51 0.941 1.269 
10 0.818 0.632 0.953 
20 0.1614 0.342 0.635 
50 0.01228 0.0832 0.262 
100 0.001577 1078 0.01623 1072 0.08382 1072 


The six transfer functions occurring in Eq. (9.2) can 
readily be calculated from Eqs. (6), (7), and (8) of 
reference 7. The result of this calculation is the 


following: 
Gia = + (2k /2)] ) 
Gina = (1/2)rC(hki) 
Gay = + (9.3) 
= (7, 4) [C(Ri) (tk, 4) | | 
Gme, = — (1/2)] 
where C(k;) is the Theodorsen function. With the 


above equations, the C. and C,, vectors have been cal- 
culated (for unit amplitude of w,(0, 0)/mo), and are 
plotted in Figs. 19 and 20. It can be seen that the 
errors in phase angle are small up to the limiting value 
of ki} = w/2 (i.e., X\/c = 2); and that the maximum 
amplitude errors, at the limiting k;, are 23 per cent 
for C, and 11 per cent for C,,.. The error decreases 


TABLE 3 
Spectrum Function ¢;(Q) for L = 1,000, «7 = 10 
LQ, La,’ = 10 = 5) = 100 
0.01 2.15 X 1077 5.93 X 107% 2.38 X 10-5 
0.02 2.15 5.93 2.38 
0.05 2.15 5.93 2.38 
0.1 2.15 5.93 2.38 
0.2 2.15 5.92 2.38 
0.5 2.13 5.92 2.38 
1 2.06 5.91 2.38 
2 1.846 5.85 2.37 
5 1.156 5.56 2.32 
10 0.482 & 107 4.89 2.23 
20 9.2). xX 3.48 1.958 
50 1.097 1.151 
100 9.44 K 107" 0.255 X 10-6 0.440 K 10% 


TABLE 4 
Spectrum Function ¢,(Q,) for L = 1,000, ¢? = 10 
LQ LQ,’ = 10 = 50) = 100 
0.01 1.1388 KX 107"! 7.45 107° 1.193 X 1077 
0.02 1.138 7.45 1.193 
0.05 1.188 7.45 1.198 
0.1 1.138 7.45 1.198 
0.2 1.138 7.45 1.198 
0.5 1.127 7.44 1.198 
1 1.108 7.43 1.193 
2 1.025 X 10-"' 7.41 1.191 
5 6.86 X 107" 7.24 1.183 
10 2.89 X 10-™ 6.68 1.159 
20 6.48 XK 107 5.17 1.072 


50 5.41 K 107 1.790 0.707 
100 6.18 X 107% 0.495 107° 0.294 107 
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rapidly as k; is reduced. As has already been remarked 


IMAG 
in Section (2), the effects of this error on the airplane 
04, response may be expected to be relatively small, since 
tae EXACT (SEAR'S FUNCTION) the gust intensity is very weak at these high wave 


numbers. 
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Laminar Heat Transfer Around Blunt Bodies 
in Dissociated Air’ 


NELSON H. KEMP,* PETER H. ROSE,* ann RALPH W. DETRA* 
lyco-Kverett Research Laboratory 


SUMMARY 


A method of predicting laminar heat-transfer rates to blunt, 
highly cooled bodies with constant wall temperature in dissociated 
air flow is developed. Attention is restricted to the case of axi- 
symmetric bodies at zero incidence, although two-dimensional 
bodies could be treated the same way. The method is based on 
the use of the “local similarity’’ concept and an extension of the 
ideas used by Fay and Riddell.'. A simple formula is given for 
predicting the ratio of local heat-transfer rate to stagnation-point 
rate. It depends on wall conditions and pressure distribution, 
but not on the thermodynamic or transport properties of the hot 
external flow, except at the stagnation point. 

Experimental heat-transfer rates obtained with correct stag- 
nation-point simulation and high wall cooling in shock tubes are 
also presented and compared with the theoretical predictions. 
On the whole, the agreement is good, although in regions of rap- 
idly varying pressure there is evidence that the local similarity 
assumption breaks down, and the theory underestimates the 
actual heat-transfer rate by up to 25 per cent. 


SYMBOLS 

¢ = mass fraction of 7th component 

= Xc(dh;/dT) 

= diffusion coefficient 

D? = thermal diffusion coefficient 

f = Eg.(7) 

= Eg. (8) 

h; = enthalpy per unit mass of ith component 

h = enthalpy per unit mass of the mixture, including dis- 
sociation energy, — h;°) 

h = heat evolved in the formation of component / at 0°K., 
per unit mass 

hp = average atomic dissociation energy times atom mass 


fraction in external flow 


H =h-+(1/2)u? 


k = thermal conductivity 

= Puktw 

L; = Lewis Number = Dipo/u 

L;? = thermal Lewis Number D;"p @,/k = 

I. = Lewis Number for atom molecule mixture 

MJ, = Mach Number of moving shock in shock tube, referred 
to speed of sound in quiescent gas in front of it 

p = pressure 

p; = initial pressure in shock tube (measured in cm. of mer- 
cury 

q = heat-transfer rate 

Ji = vector diffusion velocity, Eq. (1) 

ry = cylindrical radius of body, recovery factor 

R = body-nose radius in meridian plane 
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s = Eq. (8) 
7 = absolute temperature 
“u = x-component of velocity 
v = ycomponent of velocity 
x = distance along meridian profile of body 
4 = distance normal to body surface 
n = Eq. (6) 
0 = Eq. (8) 
uw = absolute viscosity 
= Eq. (5) 
p = mass density 
= Prandtl Number ¢,u/k 
Subscripts 
i = ith component 
e = external flow conditions 
w = wall conditions, » = 0 or y = 0 
s = stagnation-point conditions 
t,n, x, ¥ = differentiation with respect to indicated variable 
oo = free-stream conditions 
0 = local similarity solution 
E = equilibrium 


(1) INTRODUCTION 


tay THEORY of heat transfer at a stagnation point 
in a dissociating gas has been discussed by Fay 
and Riddell.' They considered the mechanism of heat 
transfer, with the pertinent physical and chemical 
effects which arise because of gas dissociation. It was 
shown that a similarity variable could be found, such 
that the equations for the stagnation-point boundary 
layer with arbitrary recombination rate could be reduced 
to ordinary differential equations. Fay and Riddell 
pointed out that the stagnation point appeared to be 
the only case in which the boundary-layer equations 
admitted this great simplification without further 
approximation. Even for the cone and flat plate, this 
reduction is only possible for the extreme cases of very 
fast or very slow atom recombination rate, correspond- 
ing to a boundary layer, either ‘frozen’ or in thermo- 
dynamic equilibrium. 

Once a stagnation-point theory is developed, the next 
step is to extend the theory to regions away from the 
stagnation point. There are certainly regimes of flight, 
such as high altitudes, where the boundary layer will 
remain laminar for some distance away from the stag- 
nation point, and a laminar theory is, therefore, of 
interest. 

While there is a great deal of literature on compres- 
sible laminar boundary layer, very little of it is appli- 
cable to the problem at hand. The distinguishing fea- 
tures of the present problem, caused by the very high 
flight velocity, are the dissociation and the large ratio 
of external to wall enthalpy (or temperature). These 


= 
| 
| 
| 
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conditions dictate large variations of fluid properties 
across the boundary layer. The variable property 
which distinguishes this from a low-speed case is the 
product of density, p, and viscosity 4, which may vary 
by a factor of 5 or more instead of being nearly con- 
stant. Most previous investigations, such as that of 
Cohen and Reshotko? and that of Levy,’ treat this prod- 
uct as a constant. Consequently, their work cannot 
be used in the present case without a careful examin- 
ation of its validity. Heat transfer by atom diffusion is 
also a new phenomenon arising in the hypersonic case. 
Romig and Dore‘ have considered a dissociating gas 
in thermodynamic equilibrium, but only for a flat 
plate with constant external flow properties and the 
particular case of equal heat transfer by thermal con- 
duction and atomic diffusion. We wish here to allow 
a variation of external properties suitable to the nose 
regions of blunt bodies, and to consider more general 
relations between thermal conduction and atomic 
diffusion. 


Lees’ has dealt with heat transfer to highly cooled 
bodies in dissociating flow. He made use of an assump- 
tion of ‘‘local similarity." At each point on the body, 
the boundary layer is assumed to be described by ordi- 
nary differential equations involving one independent 
similarity variable, with boundary conditions and 
parameters depending on local external and wall con- 
ditions. This assumption is tantamount to dropping 
certain terms in the complete differential equations, as 
will be shown below. Lees made a further approxi- 
mation by not actually solving the ordinary differential 
equations resulting from the local similarity assump- 
tion. He gave a careful discussion of the case of high 
wall cooling and concluded that the pressure gradient 
did not have a significant effect on heat transfer. He 
also concluded that the variation of pu across the bound- 
ary layer could be approximated sufficiently well by 
taking it constant at the external value. He then 
made use of the Cohen and Reshotko (pu constant, no 
dissipation) value of the nondimensional wall enthalpy 
gradient for the flat plate, thus taking it to be a con- 
stant over the whole body. 


Probstein® has proposed an extension of Lees’ work 
in which the wall enthalpy gradient is found by solving 
the energy differential equation by iteration, using 
Cohen's and Reshotko’s velocity profiles and starting 
with their enthalpy profile for pu constant. He found 
that only one or two iterations are necessary for good 
convergence in most practical cases. 


In the present work, the local similarity differential 
equations are solved exactly, with variable external 
pressure gradient parameter, variable pu, and inclusion 
of the dissipation term. Lewis Number effects are also 
discussed. This is accomplished numerically by an ex- 
tension of the method which Fay and Riddell! used at 
the stagnation point. The solutions, which depend on 
the local external flow conditions, vary around the body, 
in contrast to the constant values used by Lees. Atten- 
tion is restricted to bodies of revolution at zero angle 


of incidence. Results for two-dimensional bodies 
could be obtained by the same method. 

There are several ways to check the local similarity 
assumption theoretically. One would be to solve the 
complete partial differential equations. Another would 
be to compute the terms neglected in these equations 
in any particular case, and compare them with the 
terms retained. A third is to make use of the momen- 
tum and energy integral equations which are often used 
to construct approximate methods for predicting bound- 
ary-layer characteristics. These equations represent 
conservation of momentum and energy on the average 
across the boundary layer, and contain the nonsimilar 
terms which are dropped in the local similarity solution 
described here. In Lees’ approximate solution of the 
similarity differential equations, he takes the wall en- 
thalpy gradient constant around the body, so the non- 
similar term in the energy integral equation is iden- 
tically zero. The exact solutions of these differential 
equations presented in the present work, on the other 
hand, lead to boundary-layer characteristics which 
vary around the body, and thus afford the possibility 
of computing the nonsimilar terms in the integral equa- 
tions, and thus estimating the validity of the local 
similarity assumption in any particular case. 

In addition to theoretical checks, comparison with 
experimental results is presented in this paper for both 
hemisphere cylinders and flat-nosed bodies. The ex- 
perimental results were obtained in shock tubes at high 
stagnation enthalpies. The technique employed was 
the same as used for the stagnation point measurements 
discussed in reference 8. Only a brief description of 
the experimental techniques is given. 


(2) BouNDARY-LAYER EQUATIONS 


The boundary-layer equations suitable for the axi- 
symmetric flow of dissociating air over a body of revolu- 
tion are given in reference 1. The air is represented 
as a mixture of “air molecules’’ and ‘‘air atoms,’’ the 
differences between oxygen and nitrogen being ac- 
counted for by using average properties for the air 
particles. The diffusion velocity, measured with refer- 
ence to the mass average velocity, is then taken as 


qi = —(Dj/c;) grad c; — (D;"/T) grad T (1) 


The first term is the concentration diffusion in terms of 
the mass fraction c,; of the ith component. The second 
term is the thermal diffusion (pressure diffusion is 
neglected) and the D’s are the diffusion coefficients. 

It is assumed that each component of the air is a 
perfect gas with enthalpy /;, so that the enthalpy of 
the mixture is 


h = — 


where /;° is the heat evolved in the formation of the 
ith component at 0°K. 

The usual boundary-layer coordinate system is intro- 
duced with x measured along the body surface from the 
nose and y normal to the surface, and r denoting the 
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radius of the body cross section in a plane perpendicular 
to the body axis. 

After the usual boundary-layer simplifications are 
made, the equations of mass, momentum, and energy 
become Eqs. (13), (17), and (21) of reference 1: 


(pur), + (pvr), = 0 (2) 
puu, + pu, = —p, + (pty), (3) 


pull, + = i(k Ep) H,} y + 2) 
k Cy) , + | >-(Dip k Cp) (h; 
pci/T) (hy — (4) 


Here, /J = h + u?/2, the total enthalpy, and ¢, = 
dei(dh, dT’), which is the weighted sum of the com- 
ponent specific heats, not the specific heat of the mix- 
ture. 

Eqs. (2), (3), and (4) are suitable for the boundary 
layer in thermal equilibrium, where the composition of 
the mixture is a known function of the usual thermo- 
dynamic variables (say, p and). For nonequilibrium 
situations, an additional equation for c; must be con- 
sidered [Eq. (14) of reference 1]. In the present paper, 
attention will be focused on the equilibrium boundary 
layer. Some remarks about the possible applicability 
of the results to other cases will be made later. 

We now proceed, as in reference 1, to introduce a 
transformation of the independent variables, which 
includes the Howarth and Mangler transformations, 
as suggested by Lees:° 


x 
E(x) = dé/dx = (5) 
0 


n(x, vy) = pdy, On/Oy = pur 2¢ (6) 
0 


We also choose dimensionless dependent variabies 
based on conditions at the edge of the boundary layer, 


as follows: 
u/ue=f, f= (7) 
0 


g=H/H, 6=T7/T, 


Ss = (8) 


The functions /, g, 6, and s; depend on both é and 7. 
The quantities at the edge of the boundary layer are, 
of course, functions of x only, except for //,, which is a 
constant, and is the total enthalpy of the inviscid flow. 
(We assume that the flow at the edge of the boundary 
layer came through the normal shock near the axis of 
revolution.) 

With these transformations, the continuity Eq. (2) is 


po = + V 2€f,| (9) 


In addition, we obtain the following transformed dif- 
ferential equations, with Prandtl Number o = ¢,y/k, 
Lewis Number L; = Djpé,/k and 1 = pp/ pwbu- 


Momentum: 


+ Sfon + 2(d In u,/d In —) X 
[(p./p) = = 2E( (10) 


Energy: 
— 1)s;, + Li's + X 
— o Son} = (11) 
Eqs. (10) and (11) are the system to be solved for / and 
g as functions of £ and 9, for an equilibrium boundary 


layer. Boundary conditions are: 
7=0: f=f'=0, g= zg, = h,/H. (12) 
nao: gol 


True similarity solutions of these equations require 
all terms to be independent of £, so that they reduce to 
ordinary differential equations in n. One way that this 
can occur is for the external flow and wall enthalpy to 
be independent of &—i.e., of x—-which is true on a cone 
of constant surface temperature. 
case for high-speed flow is that of an axisymmetric 
stagnation point, at which Eqs. (10) and (11) also be- 
come similar. The external velocity u, may be written 
4, = x(du,/dx), and from Eq. (5) we easily find that 


A more interesting 


~ 


B = 2(d1n u,/d In = 1/2 (13) 
In this case, Eqs. (9) and (10) reduce to 
+ San + (1/2) — f,?] = (14) 


+ + — hi) X 
(Li — 1)s,, + Lis 6,/0|}, =O (15) 


These equations (with L;’ = 0) are the ones that Fay 
and Riddell have solved for the equilibrium case in 
reference 1. 

The purpose of the present work is to derive informa- 
tion about laminar heat-transfer rates on blunt bodies 
away from the stagnation point by reducing Eqs. (10) 
and (11) to ordinary differential equations, like Eqs. 
(14) and (15), and then solving them in the same man- 
ner as the latter equations were solved in reference 1. 
This is accomplished by the introduction of the idea of 
“local similarity.” 


(3) LocaL SIMILARITY 


When boundary-layer characteristics, especially heat- 
transfer rates, are desired for bodies with varying ex- 
ternal flow properties, we are confronted with the full 
Eqs. (10) and (11). A reasonable approximation which 
reduces them to ordinary differential equations is 
clearly extremely useful, since partial differential equa- 
tions are difficult to handle numerically. Such an 
approximation is that of local similarity, as discussed 
by Lees,® and by Fay and Riddell.' In this approx- 
imation, at any point x, the dependence of the de- 
pendent variables on £ is taken to be such that their 
derivatives, with respect to &, may be neglected. There- 
fore, the right-hand sides of Eqs. (10) and (11) are taken 
to be zero. Further, the terms on the left which depend 
on &, arising from the external flow or wall conditions, 
are assumed to have their local values. Then, the equa- 
tions again become ordinary differential equations in 7 
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Fic. 1. Heat-transfer parameter correlation. 


at any point on the body, with parameters depending 
on the local external and wall conditions, and can be 
solved for the boundary-layer characteristics at that 
point. 

Essentially then, local similarity represents a patch- 
ing together of local solutions; the x-wise history of the 
flow is ignored, except as it is contained in the external 
and wall conditions, which form the coefficients of the 
differential equations. The validity of this approxi- 
mation depends ultimately on the fact that the external 
flow properties vary slowly with é, and the terms neg- 
lected in the differential equation are really negligible 
compared to those retained. One way to determine 
this for any particular case would be to actually com- 
pute the right-side terms in Eqs. (10) and (11), once a 
solution neglecting them had been obtained. 

Another way to derive some theoretical information 
about local similarity is to integrate Eqs. (10) and (11) 
on 7 across the boundary layer from » = 0 ton = &, 
thus obtaining the boundary-layer integral equations. 
If this is done, the results are 


-{ — f,)dn + 2(d In u,/d In X 
0 


[(p./p) — + 2&(d ae) f,1 — f,)dn (16) 
0 


[(Li — + = 


— gidn + 2&(d, a) f,l — g)dn (17) 
0 0 


The nonsimilarity term in each of these equations is the 
last one on the right, which could be computed in any 
given case once a similarity solution had been obtained. 
Qualitatively, if this term is found to be small, the 
approximation of similarity certainly is a good one, 
while if it is large, the approximation has broken down. 
However, a quantitative theoretical criterion for the 
size of the term is not known. The ultimate practical 
test is experiment, and comparison of theory with ex- 
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periment should help provide such a quantitative cri- 
terion. 

Notice that if the local similarity solution does not 
depend on the external flow properties, the nonsimilar 
terms in Eqs. (16) and (17) are automatically zero, 
and there is no possibility of evaluating the error in- 
volved in the similarity assumption. Such is the case 
with Lees’ value of g,,. in reference 5. 

Once the local similarity approximation is accepted as 
one which may give useful results, the equations to be 
solved are Eqs. (10) and (11) with zero right-hand sides. 
We will neglect thermal diffusion (L;" = 0), and, for the 
time being, concentrate on the Lewis Number unity 
case (L; = 1). Actually, the Lewis Number is near 
1.4 according to best present estimates, and a way to 
take this into account will be mentioned later. With 
these restrictions and recalling Eq. (13), the momentum 
and energy equations become 


+ + — f,?] = 0 (18) 
(lg,/0), + + (u.2/H,) font, =O (19) 


with boundary conditions (12). 


(4) SOLUTIONS OF THE EQUATIONS 


The solution of the local similarity Eqs. (18) and 
(19) is based upon the solution of the stagnation point 
Eqs. (14) and (15) (with Z,” = 0) presented in reference 
1. They both satisfy the same boundary conditions, 
Eqs. (12). 

There are three points of difference between the 
stagnation point and local similarity equations. The 
first is the fact that the pressure gradient parameter, 8, 
is no longer 0.5 when one moves away from the stag- 
nation point. This difference is easily remedied by 
solution of Eqs. (14) and (15) with 1/2 replaced by other 
values. Such solutions for 8 = 0, 0.5, 1.0, and 2.0 are 
shown as the open points in Fig. 1. This plot is the 
same as that of Fig. 2 of reference 1. There it was 
found that if V 26 (1 — gy) was plotted against the 
pu ratio across the boundary layer, the re- 
sults fell on one curve, regardless of g, and /7,. This is 
seen to be true for all values of 6 in Fig. 1. In per- 
forming these calculations, the same expressions for / 
and p,/p were used as described in the Appendix of refer- 
ence |, namely 


l= pu ‘Publw = iV g — (20a) 
pe/p = 1 — — g) — y2(1 — g)4 (20b) 


The justification for using these relations will be found 
in reference 1, especially in Fig. 1 there. 

For purposes of computation, straight lines have been 
fitted to the points of Fig. 1. They are approximately 
valid in the range 0.15 < pette/ pute < 0.55, which is 
the region of interest for highly cooled walls. This fit 
yields the relation 


V — gy) = 0.648(1 + 0.096 VB) 
(Pelte/ Pwhtw) 


0.438 (2 1 ) 
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and shows the dependence on 6 of the enthalpy gradient. 
Note that the curve for 8 = 1/2 differs slightly from 
that given in reference |, but the difference is negligible 
in the range of applicability of Eq. (21). 

Having found the pressure gradient effect on the 
stagnation-point equations, let us now go on to the 
other two added effects which appear in the local simi- 
larity equations. The first is, of course, the presence 
of the dissipation term which appears as the last one 
in the energy Eq. (19). This term, which introduces 
the dissipation parameter u,”/J/,, is zero at the stagna- 
tion point. 

The second added effect is a more subtle one, involv- 
ing the fluid properties, as given in Eq. (20). There 
they are taken to depend on g, the dimensionless stag- 
nation enthalpy. This is clearly correct at the stagna- 
tion point. However, elsewhere they should depend 
instead on the dimensionless static enthalpy /, related 
tog by 


h/H, = g — f,°u.?/2H, (22) 


Thus, 4/77, must be substituted for gin Eq. (20). The re- 
sulting expression for / must also be corrected, of course, 
to have the value unity at the edge of the boundary 
layer, h./H, = 1 — u,?/2H,. Notice that this modi- 
fication of the fluid properties also depends on the dis- 
sipation parameter u,*///,. The fact that both the 
dissipation term and the fluid property modification 
depend on the same parameter greatly simplifies the 
numerical calculations. 

The size of this dissipation parameter u,°///, is 
limited. For expansion to zero temperature it ap- 
proaches 2. On practical bodies in the shock tube it 
may reach 0.5 to 0.6, and in flight the extreme value 
which could be expected is not much more than unity. 

The results of calculating some sample solutions of 
Eqs. (18) and (19) with the fluid property and dissipa- 
tion term modifications are shown in Fig. 1. First the 
fluid property change alone was inserted, for u,?///, = 
1/4, 1/2, 1.5 at several points in the 6 = 0, 1/2 lines. 
These results are shown as the /i//ed points in Fig. 1. 
Each filled point is obtained from a stagnation point 
solution modified by replacing g in Eqs. (20), by h, H, 
from Eq. (22). It is easy to see that this replacement 
increases the value of p,u,./ pwkw, SO the point is moved 
to the right. It turns out also to increase the value 
Of £4. SO the point is also moved upwards. Both these 
movements increase as u,”///, increases. The stagna- 
tion point solution which is modified in each case has 
an arrow on it, and the filled points which move up- 
ward and to the right from the arrow represent in- 
creasing values of u,*///,. 

Now to these calculations the effect of the dissipation 
term is added (for o = 0.71). The results of this addi- 
tion are plotted as the half-filled points in Fig. 1. They 
represent complete solutions to the local similarity 
Eqs. (18) and (19). It is seen that the effect of adding 
the dissipation term is to move the points straight down 
at a fixed py ratio. The important thing to note 
is that for 0.15 < pette/ Pyby < 0.6 the dissipation term 


moves the points down approximately the same 
amount as the fluid property modification moved 
them up. 

The net result of both modifications to the solutions 
of the stagnation equations is only a small change in 
the value of g,./(1 — g.), even for the largest value of 
u.”/H,. For u,?/H, = 1.5, the changes range from 
+1 per cent for high wall cooling (g, = 0.0164) to 
—6 per cent at lower cooling (g,, = 0.1682). It would 
thus appear, at least for values of the Prandtl Number 
near 0.71, that the value of g,,,. (1 — g,) is only slightly 
affected by the value of u,*//7,, and mainly by the pres- 
sure gradient parameter 8. Thus, at any point on the 
body, a useful approximation is to use the value of 
(1 — appropriate to the /oca/ 8 and the stagna- 
tion value of the pu ratio especially for values of u,* /7, 
below unity and for high wall cooling. This inde- 
pendence of the pu ratio represents a considerable 
simplification of the heat-transfer calculations, because 
the values of p and uw in the hot external stream only 
have to be calculated at the stagnation point and 
nowhere else. 

The effect of the dissipation term itself is of interest 
in connection with the so-called recovery factor, which 
may be defined here as 
r = 1 + [Agye/(Sqw)no diss.) — gu) /(ue?/2H.)] (23) 


where Ag,,. is the change in g,,, when the dissipation 
term is added to the energy equation. Sixteen cases 
were calculated for 8 = 0, 1/2 and u,2/H, = 1/4, 
1/2, 1.5 for o = 0.71. The value of 1 — r varied from 
0.145 to 0.166, or a maximum deviation of 7 per cent 
from a mean of 0.155. 

The low-speed form of 1 — r is usually taken as close 
tol — Vo. In the present case, | — 
It, therefore, seems quite satisfactory to take the re- 


o = 0.157. 


covery factor as V o for these highly cooled cases also. 
A check of one case, at ¢ = 0.5, showed 1 — r = 0.296 
compared to | — Vo = 0.297, which verifies this con- 
clusion. 

Up to now the Lewis Number for all species has been 
taken to be unity. Calculations of the Lewis Number 
effect were made in reference | where it was shown 
that for L ¥ 1 (the Lewis Number of all species was 
taken the same) a correction factor of the form [1 + 
(L* — 1)hp H,| should be used to multiply the L = 
| heat-transfer rate at the stagnation point. a@ was 
found to be 0.52 for equilibrium flow and 0.65 for a 
frozen boundary layer. /p is the average dissociation 
energy per unit mass of atoms times the atom mass 
fraction in the external flow. It seems reasonable to 
believe that a similar factor will hold away from the 
stagnation point, provided hyp is the local enthalpy in 
dissociation in the external flow. For a typical shock 
tube case of shock Mach Number 12, initial pressure 
1 cm. Hg., a check was made of the size of this factor 
using National Bureau of Standards data for equilib- 
rium air, reference 7. The ratio between the factor at 
the stagnation point, and its value after a Newtonian 
expansion to the original flow direction was 1.03 for 
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L = 1.4. Since the Newtonian law probably over- 
estimates the expansion, this represents an upper limit 
to the change in Lewis Number factor of 3 per cent. 
It would, therefore, appear that the value of the cor- 
rection factor at the stagnation point is a satisfactory 
one to use for engineering purposes over the whole body. 
Any more detailed study of this point, involving actual 
solution of local-similarity equations for L # 1, pre- 
sents a large computing problem which does not appear 
justified at this time. 

The question of heat transfer when the boundary 
layer is not in equilibrium can also be clarified by refer- 
ence to the stagnation-point results of reference 1. 
There it was shown that, although the mechanism was 
different, the total heat-transfer rate was substantially 
the same for equilibrium flow, frozen flow, and all 
intermediate cases. There is no reason to expect this 
result to change as one moves around the body. There- 
fore, it appears that the equilibrium results of this 
paper may be used even for nonequilibrium thermody- 
namical situations, provided the wall does not inhibit 
recombination of atoms. 


(5) CALCULATIONS OF HEAT-TRANSFER RATE 


In line with the local similarity ideas presented 
above, the method of calculating laminar heat-transfer 
rates will be described here. The heat-transfer rate 
to the wall g is given by the sum of the conductive and 
diffusive transport of heat, the latter being included 
only when atoms recombine on the wall. 


q = [R(OT/Oy)],=-0 + | Delhi — hi®) X 

+ Di"(c:/T) (OT /dy) |} ,-0 (24) 
As shown in Eq. (41) of reference 1, this becomes, in 
nondimensional variables, 


= 0) + DX — X 
[(L, 1)Sj, +- Li's, (25) 
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Near the stagnation point u, = x(du,/dx),, r = x, 
and = so that 


2E = V 2prsttus(du,/dx), (26) 


For cold walls—that is, where no dissociation occurs 
the summation term in the heat-transfer rate vanishes, 
and the use of Eq. (26) in Eq. (25) gives the expression 
for gs. 

Since it has been shown that effects as Lewis Number 
are accounted for by the stagnation point behavior, a 
convenient way of finding heat transfer is to calculate 
the heat-transfer distribution g/g,. For a wall on 
which no recombination occurs, for example, Eqs. (25) 
and (26) give 


’—(1/2) 


= (1 V 2¢ 2 Znws 


(27) 


To evaluate g/q, from a formula like Eq. (42) is quite 
simple using the results of the previous section. The 
velocity and wall conditions around the body must be 
found, as well as the distribution of pressure gradient 
parameter 6. Then the curve-fit of Fig. | can be used 
to find by use of 


— gue) — = + 0.09608) + 
(1+ 0.096V 0.5) = (1+ 0.096V 8) 1.068 (28) 


The pu term drops because, as explained above, the 
value of g,/(1 — g ») at any point is approximately 
the same as the value at the /ocal B and stagnation pp 
ratio. The stagnation-point velocity gradient must 
also be obtained from the inviscid flow. For example, 
for Newtonian flow, 


(du,./dx), = V 2(p, — pa)/p,/R 


where & is the body radius of curvature at the stagna- 
tion point in a meridian plane. Using this expression 
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Fic. 2. Pressure distribution on hemisphere cylinder. 
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Fic. 3. Heat-transfer distribution on hemisphere cylinder. 


and Eq. (28) in Eq. (27), g/g; can be found as a function 
of x. qs is obtained from Eq. (63) of reference 1, for 
equilibrium flow, or Eq. (65) for frozen flow. 


(6) RESULTS OF CALCULATIONS 


As one example of the use of formulas (27) and (28), 
a case for a hemisphere cylinder in a shock tube was 
calculated. The case chosen was a shock Mach Num- 
ber, J, = 12, and an initial pressure, p; = | cm. Hg. 
Equilibrium stagnation conditions are then p, = 2.242 
X 107 dynes per cm.*, 7,, = 6,945°K., (equivalent to a 
flight velocity of 18,000 ft./sec. at 70,000 ft. altitude). 
The ‘“‘free-stream’’ to stagnation pressure ratio is 
p./ps = 0.1087. The wall temperature was taken 
as a constant 300°K. 

The body was assumed to have a Newtonian pressure 
distribution, p/p; = 1 — (1 — p./ps) sin? (x/R). 
Fig. 2 shows an experimental pressure distribution ob- 
tained on a hemisphere-cylinder in the shock tube using 
Mach line measurements. The results from a number 
of experiments indicate that this formula fits the data 
within the limits of experimental accuracy. 

Three different thermodynamic assumptions were 
used to find the inviscid flow: (a) thermodynamic 
equilibrium corresponding to NBS data;’ (b) constant 
y to represent equilibrium, y = 1.126; and (c) con- 
stant y to represent gas composition frozen at the 
stagnation point, y = 1.38. 

The results for the heat-transfer distribution calcu- 
lation according to Eqs. (27) and (28) are shown in 
Fig. 3. Only the curve for thermodynamic equilibrium 
is shown because the curves from assumptions (b) and 
(c) would be almost indistinguishable from it. This 
indicates that for expansions of the magnitude obtained 
in the shock tube (p../p, of order 0.1) the thermody- 
namic assumption has little effect on g/q;. Another 
noteworthy point about this case is that the py ratio 
varied only from 0.22 to 0.29, and the ratio g,.0/£,ws 


went from 1.0 to a maximum of 1.03 and a minimum of 
0.94, thus varying only a few per cent around the body. 

The distribution obtained directly from Lees’ for- 
mula, Eq. (15) of reference 5, is also plotted on Fig. 3. 
As might be expected from the small variation of 
Just mentioned, it is in good agreement with 
the calculations of the present method, except near 
90°. The discrepancy there can be attributed to the 
deviation of u, from the straight line used by Lees. 

The results of this example tend to substantiate Lees’ 
argument that g,,, can be taken constant over the whole 
body, at least for expansions of this size, where 8 varied 
only between 1.12 and 0. Thus, Lees’ theory may be 
expected to give values of g/g, which are in good agree- 
ment with those predicted by Eqs. (27) and (28). 
However, it should be pointed out that if the heat- 
transfer rate q itself is desired, the stagnation point 
value g, should be taken from reference 1, not from 
Lees’ work. As shown in reference 10, Lees’ formula 
for g, is too low by a factor (p,4./ puby ).°" in addition 
to having no Lewis Number term. For the example 
of Fig. 3, (peu, Pulte)s’ | = 0.876 and the Lewis Number 
factor for L = 1.4 and thermodynamic equilibrium is 
1.078. Including a slight difference in the constants, 
the value of g, is thus larger than Lees’ by a factor 1.3. 

In order to test the effect of a more rapidly varying 
pressure distribution, the ratio g/q, was also calculated 
for the flat-nosed body shown in Fig. 4. Shock tube 
conditions corresponding to stagnation point simula- 
tion of 14,000 ft./sec. at SO,000 ft. altitude were used 
(M, = 9.5, pi = 1 em. Hg, simulated equilibrium y of 
1.195.) The pressure distribution used is shown in 
Fig. 4. It was obtained in the shock tube from Mach 
line measurements and was faired into the subsonic part 
of a wind-tunnel distribution obtained at a Mach Num- 
ber of 4.0 on a similar model. The velocity gradient 
at the stagnation point was adjusted according to the 
analysis given in reference 9 because the accuracy of 
the wind-tunnel measurements near the stagnation 
point leave it in doubt. 

The calculation was done by means of a program for 
the IBM 650 computer which finds g/g, from Eqs. 
(27) and (28), as well as all the pertinent inviscid flow 
properties for an arbitrary body shape. The pressure 
distribution and stagnation conditions, as well as the 
body shape, are the input data. The program uses 
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the empirical expression for gy/Z,ws, given in Eq. (28), 
and an effective value of y, which is also a piece of input 
data and may be chosen to suit the situation. 

The results of this calculation are shown as the solid 
line in Fig. 5. Note the heat-transfer peak predicted 
just at the beginning of the corner. It was found that 
the pressure gradient parameter, 8, which is extremely 
sensitive to the pressure distribution had two peaks, 
both with values of 8 = 2.35. The first peak occurred 
at the beginning of the corner, the second—halfway 
round, a little aft of the sonic point. 

Also shown on Fig. 5 is the heat-transfer distribution 
calculated from Lees’ theory by Eqs. (12) and (12a) of 
reference 5. It was obtained using the same pressure 
distribution and inviscid flow data as that used for the 
calculations by the present theory. It predicts the 
same pressure peak as the present theory but is some- 
what higher on the rest of the corner and the side of the 
body. Both theories are identical on the front face. 
There would, of course, again be a difference in the 
stagnation-point values which would make the absolute 
value of the heat-transfer rates differ. 

Experimental points obtained from shock tube ex- 
periments are shown in Figs. 3 and 5. It is seen that 
the agreement with theory is reasonably good, except 
near the corner on the blunt body where the experi- 
mental heat-transfer peak is 25 per cent higher than the 
theoretically predicted one. Further discussion of 
these experimental results is given below. 


(7) EXPERIMENTAL RESULTS 


The experimental techniques involved in making 
heat-transfer measurements in shock tubes have been 
described in a previous report.*!! 

Shock tube heat-transfer measurements have been 
made successfully with two types of heat-transfer 
gages, the thin resistance thermometer and the calo- 
rimeter gage. For the present experiment, it was felt 
to be important that the boundary layer be disturbed 
as little as possible by the presence of the gage element 
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in order to maintain laminar flow over the body. As 
a result, thin gages were thought to be superior to the 
calorimeter gages. However, similar measurements 
made with calorimeter gages, where the disturbances 
created by the gages themselves may have introduced 
some uncertainties, indicate agreement within the 
scatter. 

The heat-transfer distribution measurements re- 
ported here were performed by comparing heat-transfer 
rates measured by two gages during the same experi- 
ment. One gage was always at the stagnation point 
for reference. The points in Fig. 3 show data from 
such experiments compared to the theoretical line calcu- 
lated by the methods of this paper. One set of data 
covers the shock Mach Number range of 10.9 to 13.5 
(stagnation point simulation of 75,000 ft. altitude and 
flight velocities from 17,000 to 20,500 ft./sec.); the 
other group of experiments is performed at lower heat- 
transfer rates at a shock Mach Number of around 8.0 
(stagnation point simulation of 90,000 ft. at 12,000 ft. 
sec.). 

The data and the calculated heat-transfer distribu- 
tions of Fig. 3 appear to show that the average of the 
measurements lies somewhat below the predicted value. 
This trend is equally apparent for the three types of 
measurements shown: (1) the calorimeter gage; (2) 
the thin gage at high heat-transfer rates (J, = 12); 
and, (3) the thin gage at low heat-transfer rates (UV, = 
8). Three groups of data are presented because of 
certain experimental difficulties which are peculiar to, 
and different for, the three groups. These difficulties, 
although they most probably introduce some scatter, 
appear to have been properly accounted for. A de- 
tailed discussion of the problems can be found in refer- 
ences 8 and 11. 

In addition to the uncertainties introduced by gage 
signal interpretation, some scatter can be attributed to 
inaccuracy in the measurement of gage location (data 
is plotted at the average angle while the gage actually 
covers between 5° and 10° on the model), and also 
to the difficulty of matching the small hemisphere to 
the cylinder portion of the models very accurately. 
The latter effect is noticeable at body angles over 70°. 

In addition to the hemisphere cylinder experiments, 
the heat-transfer distribution on a flat-nosed body with 
a corner radius equal to 1/4 of the cylindrical radius has 
been measured. The small corner radius being a more 
stringent test of the similarity type of solution than 
was the hemisphere cylinder, considerable effort was 
expended in obtaining as much detail of the heat- 
transfer distribution over the corner as was possible. 
Each model was made with five heat-transfer gages, 
including at least one on the face of the model which 
was either at the stagnation point or very close to it. 
The corner are was covered in intervals of 15°. Each 
gage covered approximately 15° of the corner arc. 
The gages were spread radially around the body so that 
the influence of the disturbances caused by one gage is 
minimized in its effect on the others. 


These measurements were performed over a range of 
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shock Mach Numbers from 7.0-10.0 at two initial 
pressures, and the results are shown in Fig. 5. The 
heat transfer appears to reach a maximum value ap- 
proximately 50 per cent higher than the stagnation 
value in the vicinity of the sonic point, which occurs 
near the 25° point of the are. 

The measurements and theory shown in Fig. 5 agree 
quite well, except on the corner where the experimental 
results are 25 per cent higher than the theory and seems 
to indicate a breakdown of local similarity due to the 
rapid change in the flow as it expands about the sharp 


corner. 


(S) LIMITATIONS OF THE METHOD 


The principal theoretical assumption made in the 
method proposed here for calculating heat transfer is 
that of local similarity. Because certain terms were 
neglected in the differential equations,. the resulting 
solutions do not satisfy conservation of momentum 
and energy in the boundary layer, as represented, say, 
by the integral Eqs. (16) and (17). They, in fact, satisfy 
these equations without the nonsimilar terms which 
are the last ones on the right involving d/dé. Using 
the present solutions, it would be possible to construct 
a method for obtaining boundary-layer characteristics 
which did satisfy the complete integral equations along 
the lines of the Thwaites method for low-speed flow. 

It is possible with the present method to make a theo- 
retical estimate of the validity of local similarity. 
Since g,. (and f,,.) vary around the body, the non- 
similarity terms in the integral equations can be com- 
puted in each particular case and compared with the 
other terms. Note that such a check is not possible 
with Lees’ approximate solution of reference 5 because 
he assumes g,,,. constant around the body, and the non- 
similarity term in the energy integral is identically zero. 
In any actual case, of course, g,,, will vary and the non- 
similarity term will not be zero. 

To check the similarity assumption in the case of 
Fig. 3, the energy integral Eq. (17) is written 


= — g)dn + 2&(d/dé) — g)dn 
0 0 
(29) 


The similarity solution satisfies the equation obtained 
by ignoring the last term on the right: 


The size of the term neglected with similarity assumed 
is then 


go)dn > (d dé) (g,,00/ (31) 
) 


and, as long as it is small compared to g,,..0/¢, similarity 
may be considered a good approximation. Clearly, 


the more nearly constant g,,.0, the better the similarity 
solution satisfies the energy integral equation. 


A rough evaluation was made of Eq. (31) for the 
case presented in Fig. 3. The results showed that the 
nonsimilar term was less than | per cent of g,.9/o up 
to 25°, rose to 10 per cent near 60°, and then climbed 
rapidly, reaching 17 per cent at 70°. This would 
indicate that similarity was valid up to around 60° at 
least, but might be in doubt in the region of 70°—90°. 

Comparison of the local similarity theory with experi- 
ment in Fig. 3 indicates that the theory somewhat over- 
estimates the value of g/g, compared with the mean of 
the experimental data. However, the scatter is such 
that it is difficult to estimate the limit (in x/R) of 
validity of the similarity assumption. There seems to 
be some evidence that the agreement is worse near 
80°-90°. However, experimental difficulties in that 
region, as mentioned above, may cause most of the 
discrepancy. 

On the flat-nosed body of Fig. 4, the corner presents 
a somewhat stronger test of similarity than does the 
shoulder of a hemisphere cylinder. The geometric 
problems associated with the small size of the corner 
cause more scatter in the data of Fig. 5 than in that 
of Fig. 3, but it is clear that the theory underestimates 
the heat transfer near the corner, an effect which may 
well be due to a breakdown of similarity there. Else- 
where, the agreement is satisfactory, indicating that 
similarity is valid. 

In the practical calculation of heat transfer, the pres- 
ent uncertain knowledge of the external inviscid-flow 
properties for a dissociating gas also is something of a 
limitation. However, this does not affect the method 
presented here, which describes only calculation of the 
boundary-layer characteristics once the external flow 
is known. 


(9) CONCLUSIONS 


A method is presented for calculating the laminar 
heat transfer on blunt bodies of revolution in axisym- 
metric, highly-cooled, dissociating flow. The method 
is based on exact numerical solution of the boundary- 
layer differential equations which result from a local 
similarity assumption. The solutions are correlated 
on pressure gradient parameter 8 and py ratio /, as 
shown in Fig. 1. It is found that the effects on g,,, 
caused by the presence of the dissipation parameter 
u,”/ H, in the fluid properties and in the dissipation term 
nearly offset each other. Therefore, to a useful ap- 
proximation, g,,, is different from the stagnation-point 
value only because of the local value of 8, and is not 
influenced by the local py ratio. 


The results of the theory are embodied in the rela- 
tions 


/ — (1/2) 
4s = (1 V 2E) {2 / Lows) 


(1 — gu)] (1 = gws)| = 
(1 + 0.096V8 )/1.068 


and from reference 1° 
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0.760~°" (estes)? * IT,(1 Lws) x 
(du,/dx),'/? (1 + (L* — 1) (Ap,/H,)| 


qs 


where a@ is 0.52 for thermodynamic equilibrium and 
(0.63 for frozen flow. 

The independence of heat-transfer distribution from 
the local pu ratio is very convenient, because, combined 
with the definition of 6, it means that only the pressure 
and external velocity must be calculated around the 
body, not the external density or viscosity. The former 
two are comparatively easy to obtain—through a New- 
tonian approximation, for example—while the latter 
require use of thermodynamic properties and transport 
coefficients at very high temperatures. 

Calculations using this method are compared with 
shock tube experiments on a hemisphere cylinder and 
on a flat-nosed cylinder with corner radius one quarter 
of its cylindrical radius. The calculations and experi- 
ments are in reasonably good agreement, except near 
the corner of the flat-nosed cylinder where the experi- 
mental heat transfer peak is 25 per cent higher than the 
theoretical one. This discrepancy is probably due to 
the breakdown of local similarity because of the rapid 
expansion around the corner. 

Since the method gives values of g,,. which depend 
on the external pressure gradient through £, a theoreti- 
cal estimate of the size of the nonsimilar terms can be 
made, by using the boundary-layer energy integral 
equation, for example. Such an estimate was made 
for the hemisphere cylinder case and showed that at 
70° the nonsimilar term in the energy integral was about 
17 per cent of the similar terms. 
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On Stress Interaction in Fatigue and a 
Cumulative Damage Rule’ 


ALFRED M. FREUDENTHAL* ann ROBERT A. HELLER ** 


Columbia University 


SUMMARY 


The paper presents the underlying physical concepts as well 
as the results of an investigation of the interaction of repeated 
stress cycles of variable amplitude in producing progressive 
(fatigue) damage. The materials investigated were 2024 alumi- 
num alloy and 4340 aircraft steel; and the applied load spectra 
were based on simplified exponential distributions of gust loads 
and maneuver loads, as suggested by Lundberg. The fatigue 
tests were performed on rotating bending fatigue machines of 
special design, in which randomized sequences of load amplitudes, 
derived from the given spectra, were applied. 

The results of the investigation support the proposition of a 
quasi-linear cumulative damage rule, deviating significantly from 
the simple linear (so-called Miner) rule, but equally easy to apply 
in design. They also demonstrate the lack of design significance 
of the conventional endurance limit in ferrous metals. 


SYMBOLS 


b,c = empirical constants 

e = base of natural logarithms 

g,G = functions 

h = load parameter (slope of exponential spec- 
trum in semilogarithmic scale) 

H = number of gust cycles equal to or exceeding x 

Hy = total number of gust cycles in 10° flight miles 

i = subscript indicating reference to ith stress 
level 

k = constant with dimension of stress 

ED = parameters (see p. 436) 

IN) = probability of ‘survival’: probability of 
fatigue life > N 

m = subscript referring to mth stress level 

n = total number of discrete stress levels in 
spectrum 

N = fatigue life (number of cycles to failure) in 
general 

No = minimum fatigue life (minimum number of 
cycles to failure) in general 

N = mode of statistical distribution of fatigue 
lives 

No.o1 = fatigue life at (N) = 0.01 

N = fatigue life associated with S 


The above symbols for the various fatigue lives are used with 
the following subscripts and /or superscripts: 

N,, Nos, Ns, No-ors, Or Nei, Nosi, No-orsi, refer to fatigue lives 

at constant stress amplitudes directly observed and used 
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to estimate spectrum fatigue life on the basis of the usual 
linear damage rule 

Nr, Nor, Nr, Noor, refer to fatigue lives under randomized 
(spectrum) loading estimated on the basis of the usual 
linear damage rule 

N'., N'siy refer to interaction 
fatigue lives at constant stress amplitude 

refer to fatigue lives under random- 
ized (spectrum) loading directly observed or estimated on 
the basis of damage interaction rule 


pi = relative frequency ratio of cycles of stress 
amplitude S$; in spectrum 
p(x), p(s) = probability density functions 
P(x) = cumulative probability function 
x 
P(x) = p(x )dx 
Ps) = complementary probability function 
{1 — P(x)] 
r. = area of load spectrum below endurance limit 
r, ¥5 = number of striations per unit length of crys- 
tal under constant shear stress amplitude 
r;’ = number of striations per unit length of crys- 
tal under spectrum loading 
ae a = constant stress amplitude 
Sa = stress amplitude expected to appear once in 
108 flight miles 
g = upper limit of stress interaction phenomenon 
(estimated ) 
Se = conventional endurance limit 
— = endurance limit in randomized (spectrum) 


test (estimated ) 


Suse = maximum stress amplitude used in ran- 
domized (spectrum) fatigue test 

Suits = minimum stress amplitude used in random- 
ized (spectrum) fatigue test 

a. ed = stress ratio obtained by dividing respective 

As = difference between adjacent test stress ratios 
(si+ 1— Si) 

T(x), T(s) = return number 1/P*(x) 

V., Vai, Ve’, = “characteristic” values of extreme value 


distributions of fatigue lives N,, N.i, N,', 
N, i’, Nr, Nr’, atl = 1/e 
v = fatigue life at intersection of constant ampli- 
tude and variable amplitude (fictitious) 
fatigue diagrams (Fig. 3) 


Vii’, Vr, Vr’ 


v, v’ = log Vr/V, log Vr'/ V 
x,2 = hs = random variables 
a = scale parameter of extremal distribution of 
fatigue lives 
Qi = cycle ratio pj Np’ / Ng (observed ) 
7, 6, € = empirical parameters 
a 

= gamma function 

rv) = incomplete gamma function 


v = slope of log S-log V, diagram 
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STRESS AMPLITUDE 


FATIGUE LIFE Ng Ng 


Fic. 1. Real and fictitious S-N, diagrams. 


= slope of log S-log V,’ diagram 

= standard deviation of In N 

Se = ultimate tensile strength 

T, Ti = resolved shear-stress amplitude in the oper- 
ative slip plane 

resolved shear-yield-stress in the operative 


slip plane 

@, Wj = stress interaction factors at S, S; 

= N/N’, average stress interaction factor, re- 
ciprocal of sum of cycle ratios 


(1) INTRODUCTION 


be PRESENT INVESTIGATION is concerned with the 
development of a procedure for the prediction of 
the expected fatigue life under variable stress ampli- 
tudes of the type encountered by aircraft wings, based 
on the results of conventional constant amplitude fa- 
tigue tests. It uses a set of specially designed rotating 
bending fatigue testing machines in which arbitrarily 
selected sequences of stress amplitudes can be applied 
to the specimen,' as well as a set of standard rotating 
beam fatigue machines for the comparative tests at 
constant stress amplitudes. The basis for the statisti- 
cal evaluation of the results is the three-parameter 
distribution function of smallest values;? the physical 
considerations are based on the recently developed con- 
cepts of severely localized fatigue damage in “‘persist- 
ent’’ slip bands* or ‘“‘striations,’’’ as well as on the 
hypothesis that fatigue damage progressing under re- 
peated stress cycles of constant amplitude is unlikely 
to remain unaffected by intermittent stress cycles of 
higher amplitude, while it might be roughly independ- 
ent of stress cycles of lower amplitude. This hypoth- 
esis is reflected in the observations’ that the distance 
between ‘“‘striations,’’ in which the fatigue damage is 
localized, decreases with increasing stress amplitude, 
so that even a few cycles of high stress amplitude, inter- 
mittent between repeated cycles of low amplitude, will, 
by starting the development of new striations between 
the existing ones associated with the low stress ampli- 
tude, accelerate the rate of damage accumulation 
under the subsequent low stress amplitudes. Hence, 
interaction must be expected between stress cycles 
of different amplitudes producing fatigue damage 
concurrently: the larger the difference between two 
stress amplitudes the stronger the damaging effect of 
the intermittent higher on the lower amplitude. More- 
over, the most pronounced interaction effects will un- 
doubtedly be those produced by a number of high stress 
cycles so small that their direct fatigue damage (with- 
out interaction) would be insignificant; the closer the 


number of cycles of high stress amplitude to the num- 
ber producing fatigue failure at this amplitude, the less 
significant their interaction with stress cycles of low 
amplitude, because the total number of cycles to failure 
will increasingly be governed by the high stress ampli- 
tude alone. 

It is precisely because of the fact that fatigue under 
actual flight conditions, characterized by stress spectra 
with very low frequencies of occurrence of the relatively 
high stress amplitudes, can be reproduced neither by 
the conventional constant stress amplitude fatigue tests 
nor by the more recently favored two-level tests,® due 
to their complete disregard or inadequate consideration 
of the stress-interaction effects, that the present study 
has been undertaken. In this study it is assumed, 
for the sake of simplicity, that the effect of such inter- 
action between a constant stress amplitude and all 
higher stress amplitudes in a stress spectrum can be 
expressed by an “interaction factor’ for this stress 
level, indicating the reduction of the constant ampli- 
tude fatigue life N, caused by the intermittent action 
of the high stress amplitudes (in the absence of notches, 
only damaging—that is, “‘life-reducing’’—interaction 
seems possible). Studying this effect consecutively 
for the different stress amplitudes of the spectrum, a 
set of interaction factors is obtained which define re- 
duced fatigue lives at all stress amplitudes of the spec- 
trum; plotted against the stress amplitudes, these, in 
turn, define a fictitious S-N’ diagram embodying all 
stress-interaction effects arising in a given spectrum or 
spectrum type (see Fig. 1). Comparison of this dia- 
gram with the conventional constant amplitude S-V 
diagram not only provides a clear illustration of the 
stress-interaction effect in fatigue, but also a simple 
procedure for predicting fatigue life under randomly 
varying stress amplitudes encountered in flight. 


(2) THE CONCEPT OF CUMULATIVE DAMAGE 


The close interrelation between fatigue damage under 
constant stress amplitude and the accumulation of slip 
into ‘‘striations,’’ the spacing of which depends on stress 
amplitude and is established very early in the test, has 
been conclusively demonstrated by the recent work of 
Forsyth* and Thompson;* so has been the occurrence, 
within these striations, well within the first ten per 
cent of the fatigue life, of deep grooves that might be 
identified as microcracks. 

It therefore appears not unreasonable to conclude 
that the widely assumed pronounced nonlinearity with 
respect to the number of stress cycles NV, of the process 
of fatigue damage accumulation, implying the existence 
of a very extensive period of almost no damage, fol- 
lowed by a relatively short period of rapid damage 
accumulation, does not reflect the real conditions. 
It appears that the much simpler assumption of linear 
accumulation of damage provides a closer approxima- 
tion of the actual process, at least within the extended 
period before the final, rapidly accelerated propagation 
of the macroscopic crack. Hence, a uniform rate of 
damage accumulation (1/N,) per cycle of constant 
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stress amplitude producing failure after NV, cycles is in 
itself a useful concept, although the validity of the so- 
called linear rule of damage accumulation under vari- 
able stress amplitudes based on it has been largely dis- 
credited by experiment.’ Nevertheless, this rule, which 
has been proposed by Palmgren in 1924*° and restated 
much later by Miner,’ is widely used by aircraft de- 
signers. According to this rule, the fatigue life N 
under a stress spectrum consisting of stress ampli- 
tudes S;, each of which makes up a proportion p; of the 
total number NV of cycles to failure, is governed by the 
relation 


=1 or N=1/0 (2.1) 


where NV,; denotes the fatigue lives at the constant stress 
amplitudes S; and the ratio p;/N,; = a; is known as 
the ‘“‘cycle ratio.”’ 

Recent attempts to develop a more adequate cumu- 
lative damage rule have been based on the assumption 
that the shortcomings of Eq. (2.1) are essentially due 
to the assumed linearity in a; of damage accumulation. 
Introduction of nonlinear functions in a; has, however, 
not resulted in a more reliable damage accumulation 
rule for variable stress-amplitudes. 

This is hardly surprising in the light of the recently 
demonstrated mechanism of damage accumulation 
within the slip striations. Since the distance between 
the striations is determined by the stress amplitude 
very early in the test, more or less according to the 
Yamaguchi relation for slip’ 


T— = kr (2.2) 


where ry is the number of striations per unit length, 
(r — to) the excess of the resolved shear stress ampli- 
tude 7 in the operative slip plane over its yield stress 
7, and k a constant of the dimension of stress, the 
number of striations developing at a certain resolved 
shear stress amplitude 7; will be increased by the inter- 
mittent application of even a small number of cycles of 
amplitude 7;,,, > 7; where m = 1, 2... denotes the 
number of stress amplitudes exceeding r;. Such inter- 
action between 7; and 7;4,, will undoubtedly produce an 
acceleration of the damage accumulation at stress 
amplitude 7; in relation to the excess of the actual num- 
ber 7;’ of striations over the number 7; = (7; — 70)/Rk 
associated with the constant resolved shear stress ampli- 
tude 7;, provided the material shows neither strain- 
aging nor significant recovery. It is, therefore, also 
assumed that interaction effects of stress amplitudes 
Ti-m < 7; are nonexistent or negligible. Since the 
damage accumulation depends not only on the number 
of additional striations (7; — 7;) but also on their 
width, which is a function of the number of cycles of 
amplitudes 7;;,, > 7;, as well as on the additional dam- 
age produced by these amplitudes within the principal 
striations due to 7;, the total damage-accelerating 
effect at 7; of the stress amplitudes 7;,,, > 7; is neces- 
sarily a function of the extent and shape of the part 
of the load spectrum exceeding 7;, the effect of each 


stress level increasing with increasing interval (7,,,, 

— r,) and with increasing frequency ratios p;,,,. The 
total effect on the stress amplitude 7, can therefore be 
expediently expressed in terms of an overall stress 
interaction factor w;(Piim, Ti+m) > | by which the con- 
stant stress-amplitude fatigue life .V,, is reduced to 
N’,; = N,;/w;, and the linear damage rate per cycle 
increased from (1/N,;) to (1/N’,;) = wi/Ny, With 
this damage rate Eq. (2.1) can be modified for stress 


interaction!! 


= Le Nw) 


or NW’ = N,;) 
1 


Tne relation between the applied stress amplitudes 
S; producing the resolved shear stress r, and the num- 
ber of cycles N’,; to failure at each amplitude considered 
as part of a spectrum defines a fictitious S-NV’, diagram 
associated with the spectrum, on the basis of which the 
spectrum fatigue life can be predicted by the simple 
linear rule of damage accumulation (Fig. 1). Knowl- 
edge of the individual interaction factors w; = N,, N’,, 
would permit the derivation of the S-.V’, diagram from 
the conventional S-N, diagram. Alternatively, the 
ratio 

N'/N = = = (2.4) 


reflects the overall damage accelerating effects of stress 
interaction within the whole spectrum; @ is the recipro- 
cal of the sum of cycle ratios. 

It is to be expected that the individual stress inter- 
action factors w; for the stress amplitudes of a spectrum 
increase with decreasing stress amplitude S, because of 
the larger number of stress amplitudes S;,,. > S;. 
Where a well-defined endurance limit S, exists in con- 
ventional tests, as in mild steel, it can obviously not 
remain unaffected by the intermittent applications of 
stress amplitudes S; > S, of a stress spectrum extending 
to or below S,. Thus, with stress amplitudes at or 
below S, forming the lower part of a stress spectrum, 
the conventional endurance limit can only be sustained 
for a finite number of repetitions because of the damage 
initiated by the intermittent application of the stress 
amplitudes in the upper part of this spectrum and 
propagating under the stress amplitudes in the lower 
part. The concept of a well-defined endurance limit 
of the metal would therefore appear meaningless as a 
fatigue design characteristic unless the maximum stress 
amplitude of the spectrum is itself below the constant 
amplitude S,; under all other conditions of variable 
stress amplitudes, the endurance limit, when it con- 
tinues to exist, should become a function of the stress 
spectrum. This conclusion will in fact be borne out 
by the experimental results. 

Stress interaction in fatigue is likely to decrease 
within the range of stress amplitudes at which the 
characteristic mechanism of fatigue damage within the 
sharply localized striations is replaced by one of ex- 
tensive plastic deformation, crystal breakdown, and 
distortion of structure followed by gradual exhaustion 
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of ductility and fracture (alternating plasticity). 
Within this region, stress interaction factors should 
tend toward unity without reaching it; that they are 
usually significant enough to justify their consideration 
follows from the results of the high stress-amplitude 
fatigue tests producing fatigue lives NV < 10‘ cycles; 
the linear rule of damage accumulation provides a 
sufficiently close approximation of the real condition 
only if an average interaction factor w > | is considered. 

A relatively slow asymptotic approach to the linear 
rule must be expected with the gradual transition from 
the stress spectrum to the constant stress amplitude at 
its lower limit by reduction of the severity and frequency 
of occurrence of stress amplitudes above this limit, 
which is associated with extremely long fatigue lives. 

The purpose of the present investigation is to study 
the systematic variation with the characteristics of the 
applied load spectrum of the average stress interaction 
factor @, and to establish the order of magnitude of the 
individual interaction factors w; at S; for two standard 
aircraft structural materials, 2024 aluminum and 
SAE 4340 steel, one without and the other with a defi- 
nite endurance limit. If is expressed as a function 
of the spectrum characteristics, the true fatigue life 
N’ for the particular type of spectrum can be predicted 
from that computed by Eq. (2.4). The knowledge of 
the individual interaction factors w; and of their relation 
to the stress amplitudes may provide the key to a spec- 
trum independent rule of prediction of fatigue life under 
variable stress amplitudes when, at some time in the 
future, test results based on stress-amplitude spectra 
of different types will become available. 


(3) THe Loap SPECTRUM 


The problem of simulating complex service conditions 
in a testing machine has no simple solution. Even 


JOURNAL OF THE AERO/SPACE SCIENCES—JULY, 1959 


when a testing machine that permits the application of 
an arbitrary stress sequence has been successfully de- 
signed, the problem of how to design a sequence repre- 
sentative of service conditions remains. 

The load spectrum of the airplane is a function of 
the airplane itself, of its mission, and of its operational 
environment.'? However, on the basis of a great 
number of gust records accumulated in the course of 
recent years for various types of aircraft under various 
operating conditions,'* the conclusion has been tenta- 
tively reached'! that load spectra for such different 
types as fighters and transport planes can be approxi- 
mated by simple exponential distributions, the slope 
h of which, in semilogarithmic representation, is an in- 
verse measure of the statistical dispersion of the gust 
velocities x and the only parameter of the distribution, 
1/h being the standard deviation as well as the mean 
of the frequency distribution 


p(x) = (3.1) 
The cumulative distribution P(x) = 1 — P*(x) and 
P*(x) = p(x)dx = (3.2) 


represents the frequency or probability of values ex- 
ceeding x; the return number of such values is there- 
fore T(x) = 1/P*(x). 

The standardized exponential distributions of gust 
velocities developed by Lundberg! from the observed 
gust spectra are based on a flight distance of 10° miles, 
with a total number of //, between 10° and 10’ gust 
cycles; this is the return number of the maximum gust 
intensity included in the spectrum, the equation of 
which is, therefore, 


H = H,P*(x) = 
and (3.3) 
dH/dx = Hyp(x) = Hoh e~™ 


where // denotes the expected number of gust cycles 
equal to or exceeding the gust velocity x during a flight 
distance of 10° miles. Assuming, in first approxima- 
tion, that one gust cycle produces one stress cycle, the 
variable x can be expressed either in terms of stress 
amplitudes S or, more expediently, in the dimensionless 
terms of the ratio x = S/o, = s of stress amplitude to 
static ultimate strength o,. If s, denotes the stress 
ratio for the maximum stress amplitude S, with return 
number //) since //(s,) = 1, the (negative) slope of 
Eq. (3.3) in semilogarithmic representation is expressed 
by 

h = log Hy/(sy log e) = 2.303 Ho/sy (3.4) 


With 10° < ZI) < 10’ and 0.20 < s, < 0.80, where the 
lower limit refers to transport planes and the upper to 
fighters, the slopes vary over a range 17 < h < 80." 


Three load spectra based on JJ) = 10° cycles were 
selected for this investigation (see Fig. 2) with slopes 
hy = 17.3, he = 22.9, and hy = 34.3. The values of the 


slopes were more or less dictated by the maximum 
length of the punched programing tape of the testing 
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machine, and the fact that the load applied in the ma- 
chine could not be varied continuously but only over 
six discrete steps. While slopes of the order of 50 < 
h < 80, characteristic for long-range transport planes, 
could thus not be reproduced by the testing machine 
because the length of the programing tape would have 
become unmanageable, the selected slopes are still 
within a realistic range of actual service load distribu- 
tions for civilian and military airplanes. 

To ensure that fatigue lives in the tests should not 
be prohibitively long, it was necessary to select more 
or less arbitrary lower stress amplitude limits. Equal 
intervals As between s,,;, and Sma, were selected as 
As = 0.10. Because of the introduction of a minimum 
stress-ratiO Spin > 0, the variable of Eq. (3.3) must 
necessarily be x = (s — Sin) instead of x = 5s, as values 
S < Snin cannot occur. Since it is desirable to use 
Smin aS the lowest actual test stress amplitude, rather 
than a fictitious stress amplitude at which no tests are 
performed, the proportion p; of stress-amplitude ratio 
s; applied must be specified as the proportion of the 
sum of the amplitudes between s; and s; + As = si41. 
Hence, according to Eq. (3.2), 


Sixt 
h(s-s hsmi Asi hs, 
h e h( min) d s=e min (e e 
(3.5) 


(4) ANALYSIS OF DAMAGE ACCUMULATION 


Assuming that the conventional S-V, diagram and the 
fictitious interaction S-N,’ diagram can, in first ap- 
proximation, both be expressed by straight lines in 
double-logarithmic representation, valid for fatigue 
lives VN, > N and N,’ > N, where N associated with 
the stress level S denotes the limit below which stress- 
amplitude interaction can be neglected, so that for 
N, = N,’ < N the interaction factor w = 1, Eqs. (2.1) 
and (2.3) can be expressed in terms of the stress-ampli- 
tude spectrum and the parameters of the S-V, and S-N,’ 
diagrams. The wide scatter of the observed fatigue 
lives, however, makes it necessary to associate a defi- 
nite probability level with these diagrams. 

Extensive investigations of the statistical character 
of the dispersion of fatigue test results have led to the 
conclusion that the three-parametric so-called Third 
Asymptotic distribution of extreme (smallest) values 
limited by a ‘“‘minimum life’’ Ny provides the best rep- 
resentation of this dispersion. Hence, the prob- 
ability of surviving N stress amplitude cycles is de- 
fined by the survivorship function 


= exp {—[(N — No)/(V — (4.1) 


where the ‘“‘characteristic value’ V at / = 1/e is an 
expression of the central tendency of the distribution 
and a = m/o(In N)V/6 is an inverse function of the 
standard deviation of the natural logarithm of N. 
Eq. (4.1) is valid for fatigue lives N = N, under con- 
stant stress amplitudes S as well as VN = Np under 
variable stress amplitudes. 


= 

8 

Ss 

< LOG w 

” 

N=V 


Ng- NUMBER OF CYCLES (LOGARITHMIC) 


Fic. 3. Theoretical S-V, and S-1"'. diagrams 


The derivation of the stress-amplitude spectrum 
from observations is usually based on the trend or cen- 
tral tendency of those observations; it may therefore 
be assumed that the probability level of the standard- 
ized stress-amplitude spectrum Eq. (3.3) is associated 
with the mode of an unknown distribution of spectra. 
Eq. (3.3) should therefore be combined with an S-.V, 
relation at the same probability level. Since the dif- 
ference between the characteristic value V’ and the 


mode N of the distribution? 
V-—-N=(V—-—™) ( —1 @)"4] (4.2) 


is relatively small for the most frequent values of 3 < 
a < 6, and since the actual probability level of the 
spectrum is not known, it appears expedient to select 


the S-V, relations at the characteristic value V, = J, 
for / = 1/e as representative of the central tendency. 
Hence, 

(V,/V) = (S/S) (4.3) 
and (V,'/V) = (S/S)-° (4.4) 


where v and p denote the (negative) slopes of the rela- 
tions in double-logarithmic scale. (See Fig. 3.) 
The interaction factor at the probability level 1/e¢ is 
thus defined by 


w = V,/V,' = ($/S)° while = (4.5) 
w tends toward unity either when S approaches S or 
when p approaches v. In the first case, fatigue lives 
are rather short; in the second case, stresses may be 
very low and fatigue lives therefore arbitrarily long. 
Since @ will necessarily reflect the average tendency of 
w(.S), the limiting conditions on w are reflected in simi- 
lar limiting conditions on . Thus, for very short 
lives and for very long lives, @ can be expected to 
tend toward = 1.0 but never to attain it within the 
range V< 

Combining Eqs. (2.1) and (2.3), under the assump- 
tion of continuously varying stress amplitudes and 
probability level 1/e, for which p; = p(x), Ny; = VCS), 
N’,; = V,’(S) and the sums are replaced by integration 
over s, with Eqs. (3.1), (4.3), and (4.4) the following 
expressions are obtained :' 
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~h(8~8min) ds 


1 (4.6) 
and 


(V'p/ V5) ds = 1 (4.7) 
where the ratio 5 = S/o, < 1. 

Introducing the variable z = hs and the abbrevi- 
ations A = [exp (Zmin)/(A5)"] and K’ = [exp (Zmin) + 
(h5)’|, Eqs. (4.6) and (4.7) can be written in the form 


+ 1)] = K G(y) (4.8) 


where I’, is the incomplete gamma function of limit z 
and 


dz = [P maz (p + 1) = 
+ 1)] K’'G(p) (4.9) 


where the functions G(v) and G(p) depend, respectively, 
on the slopes v and p and on limits of integration s,, ;, 
and Hence, 
@ = = K'G(p)/K = 

G(p)/G(v) (4.10) 


or Ve’ = G(v)/G(p) (4.11) 


By Eq. (4.11) the true fatigue life under variable stress 
amplitudes can be predicted from the conventional 
S-N, diagram, provided p is known as a function of 
h and of the limits of the stress-amplitude spectrum. 

Since this function is the result of the complex stress 
interaction during fatigue tests, in which an exponential 
stress-amplitude spectrum is applied in such a way as to 
avoid all sequence effects (as otherwise the independent 
probability considerations would be invalid), the im- 
mediate object of the experiments was the determin- 
ation of the overall interaction factor @ according to 
Eq. (4.10), where lV’, is obtained from the statistical 
evaluation of a sufficiently large number of fatigue 
tests under randomly applied exponential spectra, 
while Vz is computed by Eq. (2.1) on the basis of the 
conventional S-N, relation at the probability level / = 
1/e. When @ is expressed as an empirical function of 
Ve obtained from experiment, the transcendental 
equation, Eq. (4.10), for p can only be solved by trial 
and error. For a given conventional S-N, relation 
at the probability level / = 1/e and a stress-amplitude 
interval specified by S,,;, and Sy; it provides an ex- 
pression for p as a function of h and Vz with the aid 
of which the fatigue life under an exponential stress- 
amplitude spectrum of arbitrary slope h can be pre- 
dicted. 

The above analysis is applicable only as long as the 
minimum stress amplitude S,,;, is above the conven- 
tional endurance limit S,. When S, > Sy, in, the inter- 
action effect within the part of the stress-amplitude 
spectrum between S, and S,,;, can no longer be analyzed 
with the aid of stress interaction factors since the con- 
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ventional S-V, relation does not exist within this 
range, while the fictitious S-),’ diagram may extend 
over part or over the whole of this stress-amplitude 
range, depending on the severity of the reduction of 
S, toward a new level S,’ < S,, due to the damaging 
stress interaction effects of the spectrum on the endur- 
ance limit. Under these conditions the lower limit of 
integration in Eq. (4.8) is Zmin = Ase = h(S,/o,), and 
in Eq. (4.9) nin = hs,’ = h(S,’/¢,) if Se’ > Snin. The 
new endurance limit S,’ can be obtained in rough 
approximation by specially designed tests, for instance, 
by gradually raising the lower limit S,,;, of an applied 
stress spectrum, and comparing the fatigue lives ob- 
tained for different locations of S,,;, < S, with respect 
to S,. The extent of fatigue damage produced by the 
part of the spectrum between S,’ and S,,;, is determined 
by comparing the fatigue lives obtained with and with- 
out this part. 

Because of this rather complex effect of the spectrum 
itself on the limits of integration of Eqs. (4.8), (4.9), 
and (4.10), an expression for p as a function of / and 
Vp alone may not be obtainable, while its dependence 
on S, may be obscure. 


(5) EXPERIMENTAL PROCEDURE AND RESULTS 


(a) Apparatus and Specimens 


The variable stress-amplitude tests were performed 
on vertical rotating beam fatigue machines with six 
load levels.'* The stress intensities are controlled by 
six potentiometers with six relays that are opened or 
closed by signals transmitted to them from a traveling 
punched tape on which a stress-amplitude sequence 
representative of a specified spectrum has been pro- 
gramed. Special arrangements were made to permit 
the application of the two highest stress amplitudes 
by time-controlled signals independently of the tape- 
controlled signals; this makes it possible to reduce the 
frequency of application of these two levels to below 
one per tape length, as required by the spectra of 
relatively high values of h and maximum interval 
Smin)- 

The speed of rotation of the random fatigue ma- 
chines was kept at 3,600 r.p.m. and the traveling speed 
of the tape was selected so as to make the shortest 
applied stress-amplitude block, controlled by the hole 
of shortest length on the tape, equal to 12 cycles. In 
order to reduce sequence effects of the loading to a 
minimum, random time series containing the six load 
levels in the proportion specified by the individual 
spectra were derived with the aid of tables of random 
numbers. Obviously, the limited length of the punched 
tape prevents complete randomization; however, by 
performing two sets of tests of 20 specimens each, with 
two different tapes based on the same spectrum but 
containing different sequences, it was shown that the 
sequence effect produced by repetition of the tape is 
negligible, as the fatigue lives of the two sets of speci- 
mens appeared to belong to the same statistical popu- 
lation.'! 
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For the constant stress-amplitude tests, three con- 
ventional horizontal ‘‘Krouse’’ rotating bending fatigue 
machines were used. Each of these machines was 
carefully correlated with one of the random load testing 
machines to ensure comparability of the test results by 
ensuring their statistical homogeneity; this correlation 
was periodically checked. 

Carefully polished round specimens of 2024 aluminum 
and SAE 4340 steel, 3 in. long, with a 5/16-in. maxi- 
mum diameter and a central |-in. long section gradually 
reduced to a 3/16-in. minimum diameter, were used 
in this investigation. 


(b) Testing Procedure 


The S-N,-] diagrams for the two materials obtained 
from sets of 20 tests’ '’ at several stress levels are 
presented in Figs. 4 and 5. The indicated straight- 
line approximations of the presented relations are used 
as theoretical S- 1’, diagrams. 

The equations of these lines based on Eq. (4.3) are: 
for aluminum ; 


= (0.735/s)*, O<s< 0.55 | 


and (5.1) 
(0.95/s)®-*, 0.55 1 


S 

I 


for steel; 

V, 
with = 10‘ cycles, and o, = 64,000 psi and 140,000 
psi, respectively. The selection of the value of V as 
the limit between fatigue and ‘‘alternating plasticity” 
is arbitrary. The existing experimental evidence is, 
however, inadequate for a sharper delimitation of the 
phenomenon of true fatigue than 10? << V < 10‘. The 
selection of 7 = 10, the upper limit of this range, is 
probably too high; while V = 10% might be a more 
realistic value, it is below the practical range of testing 
at frequencies of load repetition that would produce 
fatigue. Because of the dependence of Eqs. (4.8) to 
(4.11) on and on and Sin, randomized variable 
stress amplitude fatigue tests using six stress amplitudes 
were performed for the three distributions summarized 
in Table | and characterized by the three values of the 


(0.854/s)", 0.5385 << 5s <1 (5.2) 


TaBLe 1 
Frequencies of Occurrence p; of Stress Ratios s; 


Test Series 
2024 Aluminum Alloy = 64,000 pei 
22 822000 145600 026640 004580 - 001000 - 00018200 
23 17.3 822000 145600 026640 - 004580 00100000 00018200 
17.3 822000 145600 026640 00458000 00100000 
26 22. 900000, 090000 009000 000900 000080 00001000 
a7 22.8 900000 - 090000 009000 . 000900 00009000 00001000 
28 22.0 . 900000 090000 009000 . 00090000 00009000 
29 22.9 900000 5 090000 - 00900000 - 00090000 
30 34.3 . 966400 . 030600 . 001000 - 000030 - 00000100 . 00000003 
31 34.3 968400 030600 - 001000 00003000 00000100 
32 34.3 968400 030600 00100000 00003000 
SAE 4840 Steel Alloy = 180,000 psi 

? 17.3 822000 145600 026640 004580 001000 00018200 

8 i7.3 - 822000 145600 026640 004580 - 00100000 

17 148600. 02 - 00459000 
10 22.9 900000 . 090000 - 009000 . 000900 - 00001000 
22.9 . 900000 090000 009000 000800 00009000 
12 22.8 009000 00090000 
13 34.3 968400 030600 000001 00000003 
34.3 968400 000030 00000100 


70r 
65} 4 4 + 
* Vs - PROBABILITY OF SURVIVAL = l/e 
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Fic. 4. Real S-N-/ and fictitious S- V’, diagrams for 2024 alumi- 
num alloy. 


PROBABILITY OF SURVIVAL 


*Vs -PROBABILITY OF SURV'VAL = 


KS! 


S- STRESS, 


RANGE oF 


2 4 6 60 2 4 6 6 08 2 « 
FATIGUE LIFE, CYCLES 


Fic. 5. Real S-N-/ and fictitous S-V', diagrams for SAE 4340 
steel alloy. 


slope 4; = 17.3, hy = 22.9, and h; = 34.3, with two 
values Of Sar, and four values of s,,,, located, respec- 
tively, at the lowest, and the second, third, and fourth 
lowest stress amplitude of the spectrum. Thus, several 
systematic combinations of /# and S,,;, were used for 
each material, with 20 test repetitions for each com- 
bination to permit statistical interpretation of the test 
results. 


(c) Test Results 

The test results are presented in Tables 2 and 3 (see 
reference 17). A summary of test data previously ob- 
tained with the aid of load spectra of nearly exponential 
shape or of exponential spectra with much shorter 
return periods of the maximum aplitude'® is presented 
in Tables 4 and 5. These results, whenever relevant, 
will be used to complement the results of the tests with 
systematic variation of h and S,, ;,. 


(6) EVALUATION AND INTERPRETATION OF TEST 
RESULTS 


The replacement, in the actual tests, of the con- 
tinuous stress-amplitude spectra by six discrete stress- 
amplitude levels s; makes it necessary, for purposes of 
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TABLE 2 
Number of Cycles to Failure N’g in Thousands for 2024 Alumi- 
num Alloy Specimens Under (Randomized) Exponential Stress 


Distributions 
Slope 17.3 | 22.0 
Ho. of Stresslevels | ‘ ‘ | ‘ 


Series 
Speciaen 


10 

118.6 
118.3 
126.7 
2 138.3 
145.4 
166.2 
477.0 116.0 
0.0 663.0 


ae. 

315.9 
347.4 
260.0 
160.0 


TABLE 3 
Number of Cycles to Failure N’x in Thousands for SAE 4340 
Steel Alloy Specimens Under (Randomized) Exponential Stress 


Distributions 
Slope» | 48.3 


Specieen to. 


2 
‘ 
1 2 
12 a 
“4 s 
18 
17 
20 
1880.0 ‘70000° 2160.0 
*Eetinated Value 
TABLE 4 
Evaluation of Parameters for Tests on 2024 Aluminum Alloy 
Specimens 
te. | Type | No. of Cycles to Failure in Thousands 
1 Ay 608.0 85.0 166.6 233.0 5.00 .274 
2 325.0 70.0 109.5 138.0 5.00 .337 
3 c, 611.0 90.0 150.5 260.0 4.50 - 246 
4 A 405.0 60.0 134.1 195.0 5.09 332 
5 8 217.0 20.0 74.1 107.0 4.88 341 
6 c 408.0 60.0 119.6 167.0 4.42 . 203 
, D 3,620.0 140.0 495.5 850.0 5.09 +137 
8 a’ 790.0 70.0 134.8 210.0 4.82 +171 
9 B 275.0 25.0 103.8 190.0 5.57 377 
10 ce 1,090.0 65.0 203.4 390.0 4.50 -187 
ae 163.0 14.0 56.5 85.0 4.00 
12 8° 129.0 9.0 45.8 74.0 4.00 355 
13 ce 161.0 30.0 62.8 109.0 4.00 390 
14 AL 1,580.0 140.0 306.0 640.0 5.82 194 
18 By 866.0 92.0 180.0 445.0 5.50 208 
16 1,870.0 150.0 285.0 580.0 5.30 
17 DL, 16, 950.0 1,000.0 6,523.0 13,000.0 7.10 385 
18 ct, 489.0 70.0 132.7 192.0 4.00 -271 
19 Ay 83.1 30.0 49.7 74.0 5.90 . 508 
20 By 59.9 20.0 37.6 60.0 6.00 628 
21 Cy 103.0 32.0 $1.4 72.0 5.45 499 
22 Exp 14,120.0 1,400.0 3,760.0 9,000.0 7.00 . 266 
23 Exp 2,810.0 200.0 479.0 780.0 5.87 .191 
24 Exp 492.0 50.0 81.8 108.0 4.28 
25 Exp 138.0 38.0 53.3 69.0 5.00 386 
26 Exp 48,120.0 000.0 -13,308.0 21,000.0 7.90 277 
27 Exp 5,931.0 500.0 1,420.0 2,700.0 7.41 +239 
28 Exp 733.0 160.0 259.0 360.0 6.32 .353 
29 Exp 174.0 35.0 1.3 118.0 5.51 410 
30 Exp 15, 100.0 700.0 4,400.0 9,300.0 8.38 +291 
31 Exp 996.0 80.0 477.0 930.0 7.17 .479 
32 Exp 209.0 65.0 116.0 170.0 6,80 =. 585 


* Unreliable Results 
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test evaluation and interpretation, to replace the pre- 
viously derived equations based on continuous expo- 
nential spectra by their approximations in terms of dis- 
crete stress amplitudes, expressing the integrals in 
terms of sums. Substituting p, according to Eq. (3.5) 
into Eqs. (2.1) and (2.3) for the probability level / = 
1/e and considering Eqs. (4.3) and (4.4), the equivalent 
relations to Eqs. (4.8) and (4.9) are obtained: 


V/Ve = K > — = Kglv) (6.1) 
1 


= K'D — = K'glo) (6.2) 


with z; = hs; and the abbreviations 
K = exp (Zmin)/(h5)", = exp (Zmin)/(AS)?, 
K'/K = 


The equivalent of Eq. (4.10) for the average inter- 
action factor @ is therefore 


@ = = g(r) (6.3) 


from which g(p) and p associated with the individual 
observed values Vp’ and ratios & are obtained by a pro- 
cedure of trial and error. These values are presented 
in Tables 4 and 5. 

The complex and rather obscure relation between p 
and the characteristic variables h, Sar, Smin, Se and v 
implicit in Eq. (6.3) and expressed by Vz can now be 
empirically evaluated, in a simplified form, by con- 
structing a set of dimensionless parameters, with the 
aid of which the test results are aligned along a rela- 
tively simple functional relation. Such a relation was 


introduced in the form 


— Smin) h’ (Ve 
(6.4) 


P= const. 'Se) (Smar 


with the unknown parameters y, 6, and € to be deter- 
mined from the test results. It should be noted that 
in Eq. (6.4) the significant parameters appear both 


TABLE 5 
Evaluation of Parameters for Tests on SAE 4340 Steel Alloy 
Specimens 
Spectrua “be % Wrote p | tle 
Ko. Wo. of Cycles to Failure in Thousands |  £4.6.3 
1 a 195.6 35.0 72.6 98.0 371 
2 o1.4 28.3 44.1 -310 
3 ° 380.0 30.0 64.0 100.0 .168 
4 D 954.0 50.0 133.9 190.0 140 
5 1,363.0 80.0 279.0 555.0 
Dy 4,178.0 230.0 796.0 2,080.0 .191 
7 Exp 9,330.0 300.0 1,560.0 2,800.0 6.8 166 
8 Exp 2,054.0 140.0 267.0 470.0 5.9 -130 
9 Exp 320.4 100.0 148.0 238.0 87.5 
10 Exp 49,220.0 700.0 3,480.0 8,200.0 7.2 O71 
il Exp 4,945.0 150.0 497.0 920.0 6.5 -101 
12 Exp 494.5 90.0 235.0 470.0 7.7. 478 
13 Exp 220.0 — 70, 000. 0* 9.8 
4 Exp 22,334.0 $00.0 2, 180.0 4,700.0 8.6 
18 Exp 108.7 16.0 330.0 630.0 8.0. 468 


*Estimated Value 
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1871.8 286.0 62.1 4.6 0736.6 603.4 194.6 4.5 1678.7 173.7 
| 2086.4 307.6 62.2 42.5 0260.6 «690.8 «196.9 46.5 | 1580.3 238.1 res 
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separately and in combination with each other and 
with »y through Iz. 

Evaluation of the test results for 2024 aluminum 
indicated that of those parameters the slope h produced 
the most dominant effect, while the effect of the other 
parameters was sufficiently reflected in I’,z. The 
resulting relation 


p = 2(h* (6.5) 


is plotted in Fig. 6 in double-logarithmic scale; it re- 
produces the test results fairly well. With the aid of 
this relatively simple relation, the slope of the fictitious 
S-’’, diagram at the probability level / = 1/e is easily 
obtained for specific conditions characterized by h and 
Vp: the real fatigue life V’p can now be computed 
simply by applying the linear law to the fictitious 
S-V’’, relation of negative slope passing through V. It 
is interesting to compare the widest possible variation 
of i and Vz with the associated variation of p. Thus, 
within the range between = 50, Ve = 10° and h = 
10, Ve = 10°, the variation of p is roughly between 
p = 10 and p = 4, the first value being characteristic 
for the so-called low-level (long life) fatigue, the second 
for high-level (short life) fatigue. The order of mag- 
nitude of the interaction factors at the lowest stress 
level of the exponential spectra appear to vary roughly 
between 100 to 300, indicating very severe interaction 
at this level. For instance, for test No. 22 with S; = 
22,400 psi and p = 7 (Tables 1 and 4), Fig. 4 shows 
wo, = V,/V’,, = 2 - 10°/7 - 10® = 300. 

A similarly simple relation for SAE 4340 steel cannot 
be expected because of the effect of the endurance limit 
and of its sensitivity to the applied stress-amplitude 
spectrum, particularly to the relative position of the 
lowest stress amplitude of the spectrum. This sensi- 
tivity is illustrated by a comparison of the test results 
with different levels of s,i,. Comparing the fatigue 
lives obtained for a given slope of the spectrum and 
(a) all stress amplitudes applied, with those for the 
same slope but (b) without the lowest and (c) without 
the two lowest stress amplitudes, the proportion 
of damage produced by the two lowest stress 
amplitudes can be identified; since stress ampli- 
tudes below a real endurance limit produce no dam- 
age (by definition of the endurance limit) the fa- 
tigue life cannot be affected by leaving out stress ampli- 
tudes below the real endurance limit. Thus, after 
compensating for the number of missing cycles of stress 
amplitudes S,; or S; and S. below the endurance limit, 
by multiplication of the test results with the compen- 
sating factors 1/(1 — p;) or 1/(1 — py — pe) (see Table 
6), no significant difference, within the range of sta- 
tistical variation, should exist between the fatigue lives 
observed under the six, five, and four stress-amplitude 
tests for the same slope of the spectrum, provided the 
stress amplitudes left out in the five and four ampli- 
tude tests are in fact below the real endurance limit. 
Conversely, a significant increase in the fatigue life 
(after compensation for missing cycles) by the elimina- 
tion of a stress amplitude is conclusive evidence of the 
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Fic. 6. Slope p of the fictitious S-V’, diagrams for 2024 alumi- 
num alloy. 


fact that this stress amplitude is above the endurance 
limit or, if a conventional endurance limit is known to 
exist above this stress amplitude, of the fact that the 
endurance limit has been lowered to a level below this 
stress amplitude. 

With the conventional endurance limit ratio of SAE 
4340 steel established at s, = 0.535, the two lowest 
stress-amplitude ratios s; = 0.35 and s. = 0.45 were 
selected so as to be below this limit, while with s; = 0.55 
the third lowest level was just above it. The test 
results which are presented in Table 6 conclusively 
support the assumption that the level of the conven- 
tional endurance limit is lowered by the application of 
a spectrum with stress amplitudes below this limit. 
They indicate that relatively little (series 10-12) or no 
damage (series 7-9, 13-15) is produced by the lowest 
stress amplitude (s; = 0.35), while considerable damage 
is caused by the second lowest stress amplitude (s. = 
0.45). Hence, the true endurance limit ratio must be 
between s; and ss, a substantial reduction of the con- 
ventional endurance limit, which is thus shown to be 
meaningless in variable stress-amplitude tests. 

The fictitious S-1’’, diagram with (negative) slope p 
therefore extends at least to the lowest stress ampli- 
tude. Although the concept of the stress-interaction 
factor w, = V’,,/ V’,; loses its meaning within the range 


TABLE 6 

Compensated Fatigue Life of SAE 4340 Steel Alloy Specimens 
for Tests With and Without the Inclusion of Stress Levels Below 
the Endurance Limit S, 


No. of Compen- 
Test Stress V'r, Compensat-_ sated Life, 
Series Slope Levels Thous. of ing Factor Thous. of 


No. h Below Cycles (see p. 439) Cycles 

7 2 1,550 1.0 1,550 

8 17.3 1 267 5.62 1,500 

9 0 168 30.86 5,184 

10 2 3,480 1.0 3,480 
11 22.9 1 497 10.0 4.970 
12 0 235 100.0 23, 500 
13 2 70,000 1.0 70,000 
14 34.3 1 2,180 31.65 68 , 997 
15 0 330 1,000.00 330,000 
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| endurance limit in only six out of the nine test series ' 
or (5- and 6-level tests) while it is above this limit in the 
d 
Fale) i3e remaining three series (4-level tests) suggests that the 
vee gy results of these latter may deviate considerably from 
> $s « any function representing the former; this was borne | 
” i Y's 4 out by the test results. It appears, moreover, that the b 
28 interaction effect below the conventional endurance = 
2 | . 
a limit was not sufficiently expressed by the ratio of the 
7 4 lowest stress amplitude to the endurance limit. There- 
° i fore, the specific area of the spectrum falling below the " 
me ! endurance limit P, is introduced as an additional sig- ic 
O¢6 8 H nificant parameter, while the effect of Iz itself is neg- 
a l p 1.7 rh [(s, Sy) (6 )) Sc 
1 3 my Se oe 5 .08 6. is plotted in Fig. 7 in double-logarithmic scale together tt 
X, h (5-°-XSmax-Smin) Pe with the test results it has been designed to fit. As it 
Fic. 7. Slope p of the fictitious S-V’, diagrams for SAE 4340 expected, the results of the 4-level tests, while follow- a 
steel alloy. ing the trend, cannot be reproduced and might prob- l 
ably follow a type of equation obtained for aluminum 
y yp q of 
Se < Si < Smin because V,; > © below the conventional because of the lack of the effect of the endurance limit fo 
endurance limit ratio s,, a knowledge of p permits the in both cases. However, the three series are not sufli- th 
expedient prediction of fatigue lives on the basis of the cient to devise such a relation for steel for spectra ee 
linear damage law. The form of Eq. (6.4) was found totally above the endurance limit. os 
to represent roughly the trend of the test results on SAE The fit of Eq. (6.6) to the test results is worse than ‘ 
4340 steel as summarized in terirs of p in Table 5. for the aluminum data. This is not unexpected be- ‘a 
However, the fact that the lowest stress-amplitude ratio cause of the interaction effect on the endurance limit 
of the applied spectrum is below the conventional which will necessarily be different for spectra with a ” 
10® 
| 
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50 Nog FOR TEST NUMBER 5, 
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Fic. 8. True random fatigue life N’z as function of linear life Ne for 2024 aluminum alloy. 
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ON STRESS 


stress range completely above the conventional en- 
durance limit and spectra with a range lying partly 
below this limit. 

In Figs. 8 and 9 the average interaction factors @ 
obtained from all variable stress-amplitude tests have 
been presented by plotting log N’, versus log Np. 
Since Vp does not reflect the joint effect of all signifi- 
cant parameters, as has been pointed out in connec- 
tion with the evaluation of p, a relatively wide scatter 
of the results must be expected and has actually been 
found. 

Results at three probability levels (/ = 1.0, / = 1/e, 
and / = 0.01) have been plotted so as to separate the 
scatter inherent in the fatigue test results from the 
scatter introduced by the inadequate combination of 
the significant parameters. Obviously these probabil- 
ity levels refer to the statistical variation of the test 
results under a fixed (modal) spectrum and do not in- 
clude the effects of the expected statistical variations 
of the spectrum itself. The fact that the plotted data 
for the three probability levels overlap indicates that 
the trend of the .V’, versus Vz relation is not very signifi- 
cantly affected by the probability level, and that the 
scatter is mainly the result of the attempt to establish 
a direct relation between IV’, and Vz, disregarding all 
other parameters. 

In spite of the scatter, the trend of the l’p — Vp 
relation is fairly pronounced. Considering that for 
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Ve — V and for Vp —~ © the relation should tend 
toward the limiting condition V’p, —~ Vp, a function 
of the form 

v’ =v — ve (6.7) 
where v’ = log (V'p/V), v = log (Vz/V) has been fitted 
to the trend of the test results. For 2024 aluminum the 
parameters are b = 2.13 and c = 0.93, and for SAE 
4340 steel 6 = 2.33 and c = 0.9. A degenerate form 
of this relation with 6 = 0 and c = O represents the 


straight-line relation 


v =v-—1 or 


log (Vpr/V'r) = log & = 1 (6.8) 


with the constant interaction factor @ = 10, which can 
be considered to delimit a ‘‘safe’’ life, since it delimits 
the range of lowest test results. 

The scatter of the V’’, versus |’, plot can be elimi- 
nated for the tests performed with exponential spectra 
by introducing a suitable combination of the significant 
parameters. Thus, Fig. 10 indicates that the use of the 
relation 


@= 6.3( VR [A(Smaz Smin)/Smin (6.9) 


will permit a fair estimate of the average interaction 
factor for aluminum. No improvement over Eq. (6.7) 
could be obtained for steel, probably because of the 
complexity of the endurance limit interaction. 

It should be noted that, according to the above equa- 
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Fic. 9. True random fatigue life N’p as function of linear life Nr for SAE 4340 steel alloy. 
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Fic. 10. Characteristic random fatigue life V’p as function 
of the significant variables for 2024 aluminum alloy under (ran- 
domized ) exponential stress distributions. 


tions, p < vy and, therefore, @ > 1.0, so that V’p < Ve. 
Since all results plotted represent the trend obtained 
by 20 repetitions of each test condition, the evidence 
for @ > 1.0 appears to be quite conclusive. Hence, 
the numerous tests in which values NV’p > Np have 
been observed under certain nonrandomized load se- 
quences must reflect effects that overshadow the pure 
interaction effects, such as work-hardening or strain- 
hardening due to the specific load sequence, or residual 
stresses due to the specific geometry (notches, rivets, 
welds) of the specimen. In fact, comparing the re- 
ported results for smooth unnotched specimens with 
those for notched specimens or structural connections 
it becomes evident that only for the latter specimen 
forms, in which severe residual stresses can be per- 
manently built up, are the reported sums of cycle ratios 
significantly larger than one, and, therefore, a < 1.0. 


CONCLUSIONS 


A quasi-linear rule of cumulative damage based on 
factors of stress interaction and a fictitious S-N, dia- 
gram produced by such interaction has been proposed 
and confirmed by random fatigue tests with exponential 
spectra for smooth specimens of 2024 aluminum and 
SAE 4340 steel. The developed relations permit the 
fairly reliable prediction of fatigue lives of such speci- 
mens under exponential stress-amplitude spectra of 
arbitrary slope and stress range. 


The damaging effect of low stress amplitudes when 
mixed with infrequent high amplitudes is demonstrated. 
In steel this effect eliminates the conventional endur- 
ance limit as a significant design value under varying 
stress amplitudes. 
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On the Rotational Motion of a Body 


Re-Entering the Atmosphere 


T. B. GARBER* 
The RAND Corporation 


SUMMARY 


The exact equations of motion of a body acted upon by aero- 
dynamic and gravitational forces are formulated, using inertial 
axes fixed in a spherical, nonrotating earth. With the assump- 
tion that the spin rate of the vehicle is constant and of such a 
magnitude that Magnus forces may be ignored, expressions are 
developed for the applied forces and moments of forces. 

After considering the nature of a typical re-entry path, the 
equations of motion are linearized, and the translational and ro- 
tational modes are uncoupled. With the assumption that the 
translational equations may be independently solved, the rota- 
tional equations of motion are reduced to two coupled linear 
equations with known variable coefficients. 

A solution of the planar case is obtained, utilizing a modified 
WKBJ approximation method. This technique is extended to 
the coupled case. 

The stability of the oscillations is examined, and the conclusion 
is drawn that the oscillations are bounded for a statically stable, 
spinning vehicle if Cy is positive. 

A brief investigation is made of the forced solutions, and 
finally a comparison of the results of this paper with those of 
previous studies is made. 


SYMBOLS 


Reference Angles 


0 = orientation angles relating the vehicle’s principal axes 
y to the inertial reference axes 

/ = orientation angles relating the vehicle’s velocity vector 
a to the inertial reference axes 

= orientaticn angles relating the vehicle’s position vec- 

: tor to the inertial reference axes 
6, | 
" = perturbation angles related to 6 and ¢, respectively 
Pp 
? = fixed reference angle, related tog 
3 = angles relating the angular position of the vehicle's 
f 


velocity vector with respect to the longitudinal axis 
Ag, g = the total angle of attack, equal to Va? + 6? 


Coordinate Axes 


X 

| ~ 

)’ > = reference axes fixed with respect to the earth and con- 
PA sidered to be inertial 
x 

| 

vy > = re-entry body central principal axes 
x, 
y. ? = re-entry body axes of symmetry 


Re-Entry Body Parameters 
B = cone half-angle 


moments of inertia about central principal axes 


I 
I,} = 
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mass of vehicle, slugs 


m = 

Se = static margin of stability, ft. (negative for stable 
vehicle ) 

61, = displacement of the center of mass from the longitu- 
dinal axis 

Aref. = body reference area 


Aerodynamic Terms 


| = lift due to a 

Lg = lift due to 8 

(0C,/da), Cra = lift curve slope 

D = drag 

Cp = drag coefficient 

Ap = drag parameter including Cy and other 
constants 

1, = lift parameter including C,q and other con- 
stants 

Ku = damping moment parameter, proportional 
to Crg 

Ky = damping moment parameter, proportional 

Kun = Ky + Ky 


Miscellaneous Terms 
acceleration of gravity (assumed to be constant) 


g 

h = altitude 

k = atmosphere density exponential 

po = sea level density 

V = re-entry body speed 

w, = spin rate about the x principal axis (a constant ) 
(‘) = differentiation with respect to time 


INTRODUCTION 


ie THE RECENT PAST there have been a number of 
studies published which deal with the problems 
associated with the high-speed re-entry of a vehicle 
into the atmosphere.'~* References 1 and 2 are 
analytical examinations of the planar re-entry case, 
while reference 3 is concerned with the angular motion 
in three dimensions of a re-entry vehicle. In the latter 
case, the results were obtained with the aid of an analog 
computer. 

This study is also concerned with the three degree of 
freedom rotational motion of a body during re-entry, 
subject, however, to certain simplifying assumptions. 

Figs. 1 and 2 indicate the various axis systems which 
are pertinent to the problem, while the variables of 
interest are listed and defined in the Symbols. 

The translational equations of motion of a body in- 
fluenced by both gravitational and aerodynamic forces 
are 
—(D/m) — g cos 7 sin [7 cos ¢ cos A — 

gsinn sin Acos¢@ — sing (1) 


VA cos = (Lg/m) + g sin H sin A — 
gsinncos A (2) 
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Fic. 1. Relationship of the re-entry body with respect to the 
inertial axes. 


Fic. 2. Relationship of the velocity vector and the body axes 
with respect to the inertial axes. 

Vé = (L,/m) + g cos 7 sin H cos A sin @ + 

gsiny sin Asin @ — gcosnceosH cos¢ (3) 
where 

(d/dt) (r cos sin 7) = V cos ¢@cos A (4) 
(d/dt) (rsinn) = V cos ¢ sin A (5) 
(d/dt) (—r cosy cos 7) = —V sing (6) 


The rotational equations of motion of the vehicle in 
terms of the orientation angles are 


6+ 7 cos — R.[1 — (61/I,) — 
(6cos + y cos @ sin [(6//J.) (R: + 1)] X 
sin Yw, = (VM, cos — (M, sin (7) 


(d/dt) (y cos 0) — — — + 
(6cos + cos [(67/J,) (R: + 1)] X 
cos yw, = sin + (VM, cos (8) 
— (d/dt) (y sin 0) — (67/I,) X 
(6 cos + cos sin X 
(y cos cos — 6sin = M,/I, (9) 


where 
w, =~ — y sind (10) 
R, = (I, — 1,)/I; (11) 
L=1,+8, <I, (12) 


TORQUES ACTING ON A SLOWLY ROTATING Bopy* 


In Egs. (7)—(9), and M, are the components 
of torque applied to the x, y, 2 principal axes fixed in 


* The rates of rotation are such that Magnus forces may be 
neglected, since only values of w, on the order of 1 rad./sec. or 
less are considered. 
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the body. If the re-entry vehicle were completely 
symmetrical, the central principal axes and the axes of 
geometric symmetry would be coincident. However, 
in general, the mass distribution of the body is such that 
the center of mass is displaced by a small amount from 
the geometric longitudinal axis. Fig. 3 illustrates the 
relationship between the principal axes and the axes of 
symmetry, for the conical body under consideration. 
From Fig. 3, the components of force along x, y, 2 are 


F, = F,, + F., € (15) 
F, = F,, (14) 
F, = F,, — F,, ¢€ (15) 


Thus, the moments about x, y, z are 
M, = —(F,)a — /,e) (16) 


M, (F,)a (6/, r€) (F:) 4 
Ky, (6 cos + ¥ cos sin y) + 
Ky, (& cos + Bcosasiny) (17) 


M, = (F,)al, + Ku(y cos 6 cos — 6 sin y) + 
Ky, (8 cosacosy — asiny) (18) 


The subscript .! of Eqs. (16)—(18) indicates that only 
aerodynamic forces are involved. Damping moments 
are introduced through Ay and Ay. In writing Eqs. 
(16) and (17), it is assumed that ¢ is a very small angle, 
and that 6/, is a small displacement. 

With the assumption that a and ¢ are small angles 
such that sin a, 8 = a, 8 and cos a, 6 = 1, the aerody- 
namic forces which are developed with respect to the 
geometric axes of the vehicle are 


(F,,)a = —D (19) 

(F,,Ja = — D (asin — Bcos (20) 

(F.,)4 = D(acosy + Bsin (21) 

where D is the drag, L,, ,, is the component of the lift 

vector in the x,, y, plane and L,,., is the component 
of the lift vector in the x,, 2, plane. 


With the adoption of an exponential atmosphere, the 
lift and drag terms have the following form: 


Us (1 2) po Aret.(OCL, ys 
(asin — Bcos (22) 


(out of paper) 


(out of 
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z 
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Fic. 3. Relationship between body axes of symmetry and the 
central principal axes. 
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| | 2) pol ret (CCL Oa), zy X 
(acosy + Bsin y)V%e%" (23) 


D = (1/2) poA ret.Cp (24) 


where pp is the sea-level mass density, and k is a nega- 
tive number with a magnitude of approximately 4 

With the re-entry body perfectly symmetrical about 
x,, the two lift curve slopes, (OC,/Oa),, ,, and (OC, + 
are identical, Ky, equals Kyy,, and Ay, equals 
Ky,. 

Utilizing Eqs. (13)-(15), (19)-(21), and (22)-(24), 
the applied torques of Eqs. (7) and (S) have the follow- 
ing form: 


M, = 0 (25) 


M, cos ¥[{1 — (6//J.)] — M. sin y = 
(1 2) purl rer. V2e""} a[(OC,. Oa) + Cp |, 
Cpél. cos + + Kya (26) 


M, sin — (6//7,)]| + M, cos = 
(1/2) peel rer, V2e*"} B[(OC_/Oa) + Cp] 1, — 
Cpél. sin + Kuy cos@ + Ky6 (27) 


Again, in writing Eqs. (25)—(27), the products of small 
quantities have been dropped. 


THE LINEARIZED EQUATIONS OF MOTION 


Unless the vehicle is subjected to very large aero- 
dynamic disturbances during re-entry, the angles re- 
lated to lateral motion, n, A, and y, are always quite 
small. The range angle // is also small, on the order 
of 5° or less, for a typical re-entry path. However, the 
angle between the velocity vector and the X, V plane, 
¢, may vary from about —20° to —90°. For a stati- 
cally stable missile, 6, the angle between the longitudinal 
principal axis and the X, Y plane, will behave in a 
manner similar to @. 

Thus, it will be assumed that the angles n, A, y, and 
II are small such that the sine of any one of the angles 
is equal to the angle, and the cosine of any one of the 
angles is equal to one. In a similar manner 


¢$=ort+® (28) 
(29) 


where ® is a fixed reference angle, intermediate in mag- 
nitude between the initial and final values of ¢. The 
specific choice of @ is dependent upon the particular 
trajectory under consideration and the drag character- 
istics of the re-entry body. 

Eqs. (1)-(6) may now be written as follows: 


= —(D/m) — gsin — g(¢p + H) cos (30) 
Vop = (L,/m) — g cos ® + g(g@p + HZ) sin ® (31) 
VA cos 6 = (Lg/m) — gn (32) 

(d/dt) (rH) = V(cos © — ¢p sin ®) (33) 


(d dt) (rn) = VA cos ® (34) 


(d/dt) (r) = V (sin ® + @p cos ®) (35) 
Combining Eqs. (33)—(35) with Eqs. (30)—(32) yields 
V = —(D/m) — gsin ® — + HW) (36) 

op = (L,/mV) — (gril V?) (37) 

Acos @ = (Ly/mV) — (g Acos ®/V sin &) + 
(grn V? sin ®) (38) 
From Fig. 3 expressions may be obtained for @ and 8; 
a= Op — op + ecosy (39) 
B = ycos ® — Acos® + esin y (40) 


If it is assumed that V and the range rate, rH, are 
weak functions of a, and if the lateral motion of the 
vehicle is primarily oscillatory, then the rotational and 
translational motion may be uncoupled by differenti- 
ating Eqs. (39) and (40) and introducing the results 
into Eqs. (36)—(38). 

V = —(D/m) — gsin (41) 
6p = a+ ( 2) pol A rer. m) (OC; Oa) a — 
V2) + esin (42) 


y cos ® = B + (1/2) (OC; Oa) B — 
ecos (43) 


The introduction of Eqs. (42) and (43), in conjune- 
tion with Eqs. (25)—(27), into Eqs. (7)—(9) yields 


& + da + Ma + w, (1 — R.) (6+ cB) =f, (44) 
B+ dB + 2B + w,(R: — 1) (a + ca) = fg (45) 


where 
d= (1 2) po( A rer. mM) (OC, Oa) — Kun (46) 
c= d + K ww (47) 
0? —(1 2) ret. x 
+ Cp] + — Kye (48) 
= —(1 2) po( Al ret. I,) cos 
(d/dt) (grH1/V?) — cos — 
Ku(gril V?) + K weé sin Yu, (49) 
fa = —(1/2) sin pV%e*" — 


e sin Yw,"R, — Kye cos wo, + 
(R, l)w,(grH (50) 


SOLUTION OF THE ROTATIONAL EQUATIONS OF MOTION 


In obtaining Eqs. (44) and (45), products of a or 8 
with 6///, have been neglected as second-order quan- 
tities. As a consequence, one general requirement for 
the stability of a rigid body rotating about its x axis 
may not be deduced from Eqs. (44) and (45). This 
necessary requirement is 


@ & S$ (51) 


which may be verified by an examination of Eqs. (7) 


and 
As an example of the method of solution that has 
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Fic. 4. Planar angle of attack envelope as a function of altitude. 


been employed, let us examine the homogeneous form 
of Eq. (44) with w, equal to zero. 


+ daé+ = 0 (52) 


The coefficients of Eq. (52) are known functions of 
time, or altitude, once the solution of Eq. (41) has been 
obtained. As a consequence, an exact solution for a 
in terms of exponential functions is rather unlikely. 
However, assuming a solution of the form 


a= 
yields the following operational equation 
d+ =0 (53) 


Obviously, with constant coefficients, Eq. (53) reduces 
to the usual characteristic equation. A somewhat more 
illuminating expression may be obtained by assuming 
\ to be a complex variable. Thus, with A = o + jw 


= —(d/2) — (1/2) (@/w) (54) 


w? = Q? — (d*/4) — (d/2) + (3/4) — 
(1/2) (@/w) (55) 


a 
(=). —[16/,m (A, + Ap)/T, sin? &k?] 4+ 18 


where hh = l'sin ® 
Kuy = 
= (1/2) po(Aret./m) (0C,/Oa) 
Ap = (1 2) po(A ret. m)Cp 
ekh 


Fig. 4 indicates the variation of the angle of attack en- 
velope with altitude for three different values of the 
half cone angle, B. Newtonian aerodynamics is 
assumed to be valid over the entire re-entry path.* 

For the planar case the forced solution, due to the 
terms of Eq. (49), is 


* The assumption of constant stability derivatives is not re- 
quired by the method of solution, but is adopted as a reasonable 
approximation for the high-speed regime being considered. 


A solution of Eq. (52) is thus 


a= (e~“ [.4; cos Sw dt + As» sin Sw dt} 
(56) 


The problem remains, of course, to determine w from 
Eq. (55). Before this can be done it is necessary to fix 
the aerodynamic characteristics of the re-entry body. 
Previous studies': * indicate that for a typical re-entry 
vehicle, the predominant term of Eq. (55) is that which 
arises due to the static margin, /,.. Thus, 


w? w,? + [(3/4) (@/w)* — (1/2) (@/w)] (57) 
where 
wn” = —(1/2)po(Aret./Zz) ((OCL/Oa) + 
An iterative solution for w? may be assumed. 
= wn® + (3/4) — (1/2) (58) 


where 7 = 0, 1,2... and w, =w,. Utilizing two iter- 
ations as the solution of Eq. (58), the following approxi- 
mate solution to Eq. (56) may be obtained: 


1/2 
2 [w,” + (1/16)k?V? sin? 
cos (Z — Zy) (59) 


where the initial value of the angle of attack is a», the 
initial rate of change of a is zero, 


S wdt, and < (1/16)k?V;? sin? 


An examination of Eq. (59) reveals one important 
difference between the results of this analysis, and those 
of references 1 and 2—namely, Eq. (59) does not con- 
tain a term in the exponential due to drag which could 
cause divergence. The neglect of g in the translational 
equations of motion in references 1 and 2 resulted in 
the apparent instability of high drag bodies at low alti- 
tudes. From Eq. (59), the envelope of the oscillation is 


1/4 
L + AMN)/2k sin ©] ekh (60) 


a, = (6/:/1,) + (cos 
t 
f (w/cos? Z) (cos Z | dt (61) 


The first term of Eq. (61) is a trim angle of attack which 
occurs due to the center of mass being displaced from 
the axis of symmetry of the body. The second term 
approximately represents the angle of attack that arises 
due to the rotation of the velocity vector. 

When considering the coupled case, Eqs. (44) and 
(45), the same general approach may be utilized in ob- 
taining solutions. Thus, 


(62) 


B= (63) 
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The introduction of Eqs. (62) and (63) into Eqs. (44) 
and (45) yields 


A(X + A2 + dd + 2%) + 
Bw,(1 — R:) (A +c) = 0 (64) 


—Aw,(1 — R.) A +c) + 


B(X + + dd + 2) = (65) 


With 4 and B constants, a nontrivial solution exists 
only if the determinant of the coefficients of Eqs. (64) 
and (65) equals zero.° 


o = —}dw/[2w (1,/1,)w,|} — 
}@/[2w (J, + 
w? F (1,/1,)wpw wn? + — 
As a first approximation of w, with w, > 0, 
w = +(/,/21,)w, + (1/2) WU, /1,)? + 40,2 (69) 


The four damped natural frequencies are thus 


+ + dr + + w,"(1 R.)* (A + c)* 0) | 27,)w, + + (1,0, = — ws (70) 
(66) 
= (1w,/2/,) — V + = (71) 
Again assuming A to be a complex variable, ; ; 
The corresponding damping terms are 

d d (1/2) d/dt) + 

2 4 Wea? + (/,0,/2I,)* 2 Vo,? + (1,0,/2I,)? Vw," + (Iw, 21.)? 

d (1,/1:)w,d (1/2) (d/dt) Vw? + 

02 = = = —_ 
2 4 V + (1,0, 27,)* 2 V eon? + 27,)* V w,” 


The solutions of Eqs. (44) and (45) are 


a = eft! [.1; cos Z; + sin Z,| + 
ef [13 cos Z3 Ay sin ay (74) 


8 = ef*" [B, cos Z; + By sin Z;| + 
ef att [Bs cos Z3 a B, sin 23] 4 By (75) 


where S widt, 23 = S 


Of the eight constants of integration in Eqs. (74) and 
(75), only four are independent, since only four initial 
conditions may be specified. 

Theoretically, the forced solutions of Eqs. (44) and 
(45) may be expressed in terms of integrals, once the 
homogeneous solutions have been obtained. How- 
ever, such expressions would be very complex, and, as 
a consequence, approximate solutions are more indica- 
tive of the nature of the forced angular motion. 

The functional form of a, and 8,may be ignored, when 
considering the constants of integration of Eqs. (74) and 
(75). Thus, with the initial altitude, 4), such that all 
of the aerodynamic terms are zero, the following ex- 
pressions may be found for the .1’s and B’s: 


Ay = [Bo — B,(0)| 
Ay = [ae — &,(0)] 
As = [a — af0)| — 
Ay & [Bo — — Ae 
By = [&(0) — do] 
Bz = [By — B,(0)] 
Bs = [Bo — 8,(0)| — By 
By & (Lw,/21,)"" [a (0) — av] — Bs 


In the determination of the .1’s and B’s from Eqs. 
(74) and (75), algebraic expressions for each of the con- 


stants in question should be obtained before allowing 
terms containing e*” to become zero. 

The stability of the oscillations is determined by the 
time integrals of Eqs. (72) and (73). Two cases will 
be considered. With w, > w, > 0), 


S odt = Sf odt = 
— (1/2) fd dt — (1/2) In Vw," + ([,w,/2I,)? (76) 


With w, > w, > 0, 


S odt = — djdt 


(1/2) In + (77) 


S odt = — 


(1/2) In + (78) 


An examination of Eqs. (76)-(78) and (46)-(47) 
indicates that the magnitudes of the oscillations are 
always bounded as long as the lift curve slope is posi- 
tive, and, of course, the missile is statically stable. 

In all cases, the .1’s and 4’s are multiplied by the 
term 1/[w,? + (J,w,/2/,)?|'/* which arises from the log 
expressions in Eqs. (76)—(78). Thus, even if the damp- 
ing were zero, the maximum amplitudes of a@ or 6 never 
exceed, assuming zero initial angular rates, ay or Bo. 
This is in contrast to the planar solution given by Eq. 
(60), where, even with C,, positive, amplification may 
occur due to the term 

For the typical re-entry body, w, is much greater than 
the spin rate, w,, in a band of altitudes extending ap- 
proximately from 50,000 to 200,000 ft.* Under these 
conditions, the expression for three degree of freedom 
rotational motion which is analogous to Eq. (60) is 


* For example, at 200,000 ft., an w, of from 10 to 15 rad./sec. 
would be representative. 
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A plot of the envelope of the total angle of attack, .1, , 
as a function of altitude is shown in Fig. 5. 

An examination of both the planar and the three- 
dimensional case reveals that the envelope of the oscil- 
lations diverges when the cone half angle, B, is 60°. 
This is due to the fact that, assuming Newtonian aero- 
dynamics, the lift curve slope for a cone is negative for 
values of B greater than 45°. 

Eq. (79) represents not only the envelope of the total 
angle of attack, but, with ae = Bj = f, = fs = 0, it is 
also the expression for .1, , aS a function of time or 
altitude. Thus, for the stated conditions, .1, , de- 
creases monotonically with time as long as C,, is posi- 
tive.* 

With the velocity vector as one axis of a reference 
system, the angular position of the body longitudinal 
axis of symmetry may be expressed as follows: 


where 


= (8/a) (SO) 


The rate at which the body rotates or precesses 
about |’ may be found by differentiating Eq. (S0). 
Again considering only the homogeneous solution with 
the initial conditions a and Bo, 


-Z, = (1,w, — Wn (S1) 


where w, > w, > 0). 

Thus, the body rotates about the velocity vector, 
in a sense opposite to that of w,, and in a spiral with a 

When the general solution is considered (still with 
&, > w, > 0), the motion of the longitudinal axis with 
respect to the velocity vector is usually very com- 
plicated. Not only does the longitudinal axis precess 
about the velocity vector, but nutation, an oscillatory 
variation in the magnitude of 1, 4, may occur. The 
particular pattern of motion followed by the vehicle 
depends upon the assumed initial conditions of the 
problem. 

When altitudes of 400,000 ft. or higher are considered, 
the dynamic pressure is very small, and, thus, w, > w,. 
In the intermediate altitude range from 200,000 to 
400,000 ft., the two frequencies are roughly the same 
order of magnitude, depending of course, on the spe- 
cific re-entry body characteristics and spin rates that 
have been assumed. Even though the two frequencies 
of oscillation are damped by different amounts, as is 
indicated by Eqs. (77) and (78), at altitudes above 
200,000 ft. al] of the damping terms are extremely 
small, and, as a consequence, Eq. (79) may be utilized 
even with w, > w, without significant error. (See 
Fig. 6.) 

* In fact C, could even be slightly negative because of the 
presence of the damping term A yyy which is always positive. 
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+ AMN)/2k sin ekh 
(79) 


The rotational motion of the re-entry body is anal- 
ogous to that of a pendulous gyro which has a spring 
constant term that increases with time. Thus, the fre- 
quency of oscillation increases until the altitude at 
which the peak dynamic pressure occurs is passed. 
At the same time, the amplitude of the oscillation de- 
creases due to both the increasing dynamic pressure 
and to the loss of rotational kinetic energy through the 
normal damping process. 


THE FORCED OR PARTICULAR SOLUTIONS 


An examination of Eqs. (49) and (50) indicates 
that the forced solutions of a and 6 arise due to body 
asymmetries, and to the planar rotation of the velocity 
vector. The latter effect is very low frequency in 
nature, and, in general, may be neglected when con- 
sidering the high-frequency oscillations of the re-entry 
body about the velocity vector. However, it is inter- 
esting to note that the rotation of V in the vertical 
plane leads to a unidirectional forced response in 8. 
The same effect has been noted in the case of projec- 
tiles fired from guns.® 

With w, > w, > 0, the following approximate ex- 
pressions may be written for the forced solutions of a 
and 


ay = [Ap (AL + Ap)] (6/,/1,) cos w,t — 
COS wet — (grH/V?) + 
eK sin wit + (1/w,2) (d/dt) (grH/ V2) (82) 


By + (61./1,) sin wt — 
eR(w,/w,)? sin w,t — eK w,t + 
— 1) (w,/@,2) (grH/V2) (83) 


where w,/. 

An inspection of Eqs. (82) and (83) indicates that the 
term which is due to the displacement of the center of 
mass from the longitudinal axis of symmetry is by far 
the most important forcing function. It is this mass 
asymmetry term, which, in the case of a nonrotating 
body, leads to a steady-state trim angle of attack. 


Total transient angie of attack ratio, (Ag B/4q, B10!) 
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Altitude, h (thousands of ft) 


Fic. 5. Total transient angle of attack as a function of altitude. 
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Such an angle of attack could cause the vehicle to de- 
part markedly from the normal re-entry trajectory. 
In the case of a spinning body, however, this trim con- 
dition tends to be “‘averaged out.” 

The various forcing functions that arise due to the 
fact that the principal axes and the body axes of sym- 
metry are skewed with respect to each other may in 
general be neglected. This is true, since, for a care- 
fully designed re-entry vehicle, the angle « is always 
very small, and furthermore, all of the skew asymmetry 
terms are multiplied by 1/w,”. 

The final forced terms are those that arise due to the 
fact that the longitudinal axis of the vehicle lags slightly 
behind the velocity vector. As was indicated previ- 
ously, and as can be seen from Eq. (83), this effect leads 
to a nonoscillatory forced value of 8 which could con- 
ceivably cause the vehicle to depart from the original 
trajectory plane. However, for a typical re-entry path, 
the lateral displacement of the body at impact is small, 
unless very high spin rates are considered. 

In view of the preceding discussion, the forced motion 
of the longitudinal axis of the vehicle with respect to 
the velocity vector is, approximately, 


(Ag, [Ap/(At + Ap)] (84) 
and (pa)s = (85) 
where a= mm = 8 = () 


CONCLUDING REMARKS 


First, some comments concerning the method by 
which the results of this paper were obtained are in 
order. In essence, the procedure that has been fol- 
lowed is a modified WKBJ method.’ The transforma- 
tion that has been employed, and the subsequent iter- 
ative process that has been utilized to solve the trans- 
formed equations, are very similar to the corresponding 
steps of the WKBJ approximation technique. 

In the treatment of the coupled a and 68 equations, 
roll symmetry must exist for the solution of the secular 
equation to be valid. Without such symmetry the 
operator, A, associated with the a equation must differ 
from the operator associated with the 8 equation, and, 
thus, it would be necessary to uncouple the equations 
before solutions could be obtained. 

A comparison of the results of the planar analysis 
of this paper with those of references | and 2 indicates 
that the differences which exist are primarily due to the 
neglect, in the latter studies, of g in the formulation of 
the equations of motion. At high altitudes and speeds 
such an approximation is of little consequence, but at 
low altitudes, after the vehicle's speed has decreased, 
g is very important. If gravitational forces are not 
included, |’ in Eq. (60) will, for high drag bodies, 
approach zero rather than some finite terminal speed. 

Since the results of reference 3 were obtained with the 
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Fic. 6. Total transient angle of attack as a function of altitude. 


aid of a computer, a quantitative comparison is diffi- 
cult. However, the nature of the motion indicated by 
Eqs. (79) and (81) is in agreement, qualitatively at 
least, with the results, obtained for the same initial 
conditions, shown in Fig. 9 of reference 3. 

Subsequent to the preparation of this paper, an 
analytical study which also considers the rotational 
motion of a missile during re-entry has appeared.” 
The method of solution utilized in reference § is similar 
to that of references 1 and 2, with the added restric- 
tions of a constant velocity, and the complete neglect 
of damping terms. Furthermore, no examination was 
made of the particular or forced solutions of the equa- 
tions of motion. 

Adopting the approximations of reference 4, Eq. 
(79) of this paper reduces to the expression given by 
Eq. (30) in the former study. In addition, the results 
shown in Fig. 6 of this study are in agreement with the 
conclusions in reference 8 with respect to the influence 
of spin rate upon the angle of attack envelope. 
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An ‘Insulating’? Boundary-Layer Experiment* 


F. H. Durgin 

Naval Supersonic Laboratory, Massachusetts Institute of Technology, 
Cambridge, Mass. 

January 16, 1959 


EAT TRANSFER to surfaces with variable wall temperature has 
been the subject of several theoretical studies. It is 
known, of course, that the transfer at any point is influenced by 
the upstream history of the flow. This fact may have consider- 
able practical value if artificial cooling is needed. It has been 
pointed out by Libby and Pallone! that if an upstream section of a 
surface is cooled below recovery temperature, the uncooled sur- 
face downstream will be swept by a flow with an energy deficiency 
relative to the fully uncooled flow. The surface temperatures 
over the uncooled portion will then remain sensibly below re- 
covery temperature for a finite distance. If, for instance, 
transpiration cooling is used, this ‘insulating’? property is 
obviously helpful. 

Since quantitative confirmation of this result was not avail- 
able, an experiment was performed to obtain such data. A 
composite cone with a copper forepart and an afterpart of balsa 
with birch veneer (for good insulation) was constructed. The 
height of the cone was 21 in. and the copper part was 10 in. The 
base diameter was 6 in. The copper nose was uniformly cooled 
to about 100°F. below recovery temperature, except for the first 
3'/2 in. from the tip, and the surface temperature over the cone 
was measured by a number of thermocouples. Reference 2 
gives a detailed account of the design, etc. 

The results are presented in Fig. 1 in the form of a normalized 
temperature 


* This research was supported in full or in pact by the USAF under Con- 
tract AF18(603)-82 monitored by the Air Force Office of Scientific Research 
of the ARDC. 
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TABLE 1 
m (flat plate) m (cone) 
Laminar 0.5 1.6 
Turbulent (1/7th power law) 0.8 1.8 


where 7, = surface temperature at junction of copper and birch 
parts, 7, = adiabatic temperature, and 7,, = 7,,(x) = measured 
surface temperature. This is plotted versus (x/l — 1) where x 
and L are measured along a ray from the tip. 

The tests were run in the Massachusetts Institute of Technology 
Naval Supersonic Laboratory wind tunnel at JJ = 3.5, a = 0, 
with R = 1.18 & 108 and 4.3 X 108 per ft. for laminar and turbu- 
lent flows, respectively. For the latter, a wire trip at x = 1 in. 
caused transition at x ~ 3.5in. Adiabatic recovery temperature 
data was obtained at 0.851 + 0.002 and 0.878 + 0.002 for laminar 
and turbulent cases, respectively. Schlieren observations and 
boundary-layer impact probe readings confirmed the type of flow 
that existed. 

The theoretical curves for laminar flow in Fig. | were obtained 
using the measured 7, distribution on the copper part (0 < 
x/L < 1) so that the effect of the uncooled tip region is included. 
It can be shown that all of the results of references 3-7 (for the 
ideal case of a uniformly cooled nose followed by a perfectly 
insulated aft section) can be put in the form 


(L/x)m/1-n du 
y= (sinn (2) 


In this form m takes the values given in Table 1; the rate of 
growth of the layer cross-section area is proportional to the mth 
power of the distance from the tip. The various theories give 
different values for n as indicated in Fig. 1; is determined by 
the temperature profile in the boundary layer. 

It is evident that the persistence of cooling (the “insulating” 
property) is appreciable even one nose length behind the copper 
section. The laminar theories are well substantiated while the 
turbulent data falls between the predictions of Rubesin® and 
Seban.? A value for » of 2/13 in Eq. (2) brings the theory and 
data into quite close agreement. 
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Simultaneous Calculation of Influence 
Coefficients and Influence Loads 
for Arbitrary Structures 


Bertram Klein 

Structures Department, Convair Astronautics, 
A Division of General Dynamics Corporation, 
San Diego, Calif. 

January 14, 1959 


o- PEOPLE have expressed interest in seeing how the author 
is able to calculate influence coefficients for arbitrary struc- 
tures by the matric method advocated in references 1-3. As 
already mentioned,' additional information is obtained when one 
seeks such data by means of the method. The purpose of the 
present note is to illustrate the procedure. A particular struc- 
ture is chosen in the illustration in order to make the explanation 
more concrete. However, the principles involved remain true 
for arbitrary structures. 

Let us consider the stiffened panel shown in Fig. 1. Notice 
that there are six movable joints, three fixed joints, twelve bars, 
and four shear panels. All together, one obtains 40 equations 
expressing the state of stress and strain in the structure. The 
corresponding 40 unknowns consist of 24 axial loads, four shear 
flows, and twelve displacements, designated by the symbols, 
P, Q, and 6. The matrix describing the elastic behavior of the 
structure is shown schematically in Table 1. This matrix has 
been set up completely automatically by the IBM 704 digital 
computer in less than half a minute, using a general program de- 
veloped at Convair Astronautics. The matrix contains 126 
nonzero elements out of a maximum number of 1,600, thereby 
being 8 per cent dense. Actual values are not shown here be- 
cause of space limitations. 

The matrix has been pretriangularized by the Witte technique® 
in less than half a minute and inverted. 35 out of the 40 equa- 
tions were pretriangularized, which is 87.5 per cent. The inverse 
is shown schematically in Table 2. It can be proved that the 
indicated, square sub-matrix in the lower left-hand corner of the 
inverse is the so-called and well-known flexibility influence coefii- 
cient matrix—i.e., displacements throughout the structure due 
to unit loads applied at the movable joints. However, as seen, 


Fic. 1. Stiffened panel. 
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this portion is but a small part of the inverse. Also included as 
sub-matrices are influence coefficients due to unit shear forces 
applied along all the bars, and all the loads due to unit loads 
applied at the movable joints and those due to unit shear forces. 
These latter loads may be designated as influence loads. An- 
other set of sub-matrices represent loads and displacements due 
to unit axial strains, such as may arise as the result of thermal 
effects. Finally. there are loads and displacements due to unit 
shear strains of the four panels. Many of the sub-matrices other 
than the flexibility influence coefficient matrix may be of interest 
in many practical applications. 

It should be realized that by the use of the Witte technique one 
is able to treat efficiently much larger structures than the one 
used in the present example. The density of the matrix actually 
decreases with increasing matrix size, and the per cent of the 
number of rows pretriangularized remains large. 
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An Orbital Energy Indicator* 


Harold E. Bamford, Jr. 
Senior Group Engineer, Boeing Airplane Co., Seattle, Wash.** 
January 23, 1959 


ABSTRACT 


An analog indicator display based upon the conservation of a 
spacecraft’s orbital energy is described. The display consists of 
a single rotating and vertically translating index and three verti- 
cal scales. Altitude, velocity, and orbital energy are indicated 
simultaneously. The interpretation of the display is illustrated 
by several examples. 


DISCUSSION 


— ANALOG INDICATOR display described herein was conceived 
in the course of an analysis of the control functions which 
must be allocated between the crew of a spacecraft and the resid- 
ual spacecraft system. The more general conclusions of that 
analysis, including a consideration of the utility of such displays as 
this, will be presented in subsequent reports. The scope of the 
present report is limited to the description of the Orbital Energy 
Indicator and to the single caveat that the criticism of any indica- 
tor display is meaningful only in relation to the task to be per- 
formed at the crew station where it is installed. 

The Orbital Energy Indicator is based upon the conservation 
of a vehicle’s mechanical or orbital energy. Only work by an ex- 
ternal force applied to the vehicle can change its total mechanical 
energy. In addition to its mechanical energy, of course, the 
spacecraft has reserves of energy in other forms, but these will 
only concern us as they power the propulsion system and thereby 
augment or diminish total orbital energy. In the display of the 
Orbital Energy Indicator, the spacecraft’s mechanical energy is 
analyzed into its potential and kinetic components. 

As can be seen in Fig. 1, the indicator display consists of 
three vertical scales and a single moving index. The index slides 
up and down the energy scale and pivots so as to indicate simul- 
taneously three points, one on each scale. This is possible due to 
the design of the scales. 

Total mechanical energy is indicated on the center scale by the 
pivot point of the index. This scale is linear, and its unit or mod- 
ulus is the potential energy of the vehicle at rest on the surface of 
its primary. This energy level, it will be noticed, is represented 
by a scale value of —1.1), since it is due to the downward force of 
the vehicle’s weight. 

Altitude is indicated on the right-hand scale by one tip of the 
index. It is measured in radii from the center of the primary. 
Thus, the minimum altitude indication is 1—the surface altitude. 
This is a reciprocal scale. An altitude of 7 is represented by a 
point on the scale one rth its length from the top. The point is 
located in relation to the top of the scale rather than to the 
bottom so that great altitudes may be represented on the scale 
above lesser altitudes. 

* The development reported here had its origin in a study carried out by 
Ritchie and Associates, Inc., of Dayton, Ohio, and jointly sponsored by 
Ritchie and Associates, Inc., and the USAF. The AF sponsorship was under 
the direction of Lear, Inc., of Grand Rapids, Mich., and in support of 
Contract No. AF 33(619)-5167, under Task No. 61940, Project No. 6190, of 
the Wright Air Development Center. Permission is granted for reproduction, 
publication, use, and disposal, in whoie or in part, by and for the United 
States Government. 

** Formerly, Senior Associate, Ritchie and Associates, Inc., Dayton, 
Ohio. 
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On the left-hand scale, the other tip of the index indicates the 
spacecraft’s velocity. The modulus of this scale is the velocity 
of a hypothetical vehicle in a circular orbit with an altitude of 1. 
This minimum orbital velocity is called ‘‘circular velocity.”’ If 
the primary is the earth, circular velocity is approximately 8 km./ 
sec. The velocity scale is an exponential one—a velocity of v 
is represented on the scale by a point separated from 0.0 by a dis- 
tance proportional to v?. 

To illustrate the interpretation of the display, let us suppose 
that the vehicle is resting at the surface of the primary. A ve- 
locity of 0.0 and an altitude of 1 can be indicated simultaneously 
only if the pivot point indicates an energy of —1.0. As we have 
seen, this is what it should indicate under these conditions. 

Now let us suppose the spacecraft to be removed an infinite 
distance from the primary—far enough for its potential energy to 
vanish. If its velocity is still 0.0, its energy should now be 0.0 
also. And, indeed, that is the energy indication for infinite alti- 
tude and a velocity of 0.0. It can be seen that this is the greatest 
energy possible with a velocity of 0.0. 

Now, finally, let us suppose the spacecraft to be accelerated in- 
stantaneously from rest to the velocity represented on the scale 
by the lozenge. This velocity can be indicated simultaneously 
with an altitude of 1 only if the pivot point indicates an energy of 
0.0. Clearly, the increment of kinetic energy exactly cancels the 
vehicle’s potential energy. To generalize, whenever the kinetic 
energy of the vehicle exactly cancels its potential energy, the in- 
dication on the energy scale must be 0.0. For any altitude there 
is a velocity which will yield this indication. Such a velocity is 
called ‘‘parabolic velocity,” since the unperturbed orbit of a ve- 
hicle whose energy is 0.0 isa parabola. The parabolic velocity for 
an altitude of 1, which is marked by the lozenge, is called “‘escape 
velocity.”’ 

It would scarcely be possible in practice to realize a parabolic 
orbit. Such a path is of significance only as the tiansition be- 
tween closed, or elliptical, orbits, whose energy is negative, and 
the more energetic hyperbolic orbits. In an orbit of the latter 
form, the index would rotate clockwise during the approach to the 
primary, and thereafter counterclockwise to yield a monotonic 
steepening of its slope. In an elliptical orbit, on the other hand, 
the slope of the index would fluctuate periodically. As the ve- 
hicle approaches perigee (assuming the primary to be the earth), 
the velocity indication would approach a maximum. Then, as 
the altitude passes through its minimum, the velocity indication 
would fall off again, reaching its minimum at apogee. Then a 
new cycle would begin. In the special case of a circular orbit, of 
course, the index would not rock at all. 

In the preceding discussion it has been assumed that total orbi- 
tal energy is constant and, hence, that the pivot point of the in- 
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dex in the display does not rise or fall. As we have noted, orbital 
energy can only be changed through the action of an external 
force. During the spacecraft’s initial ascent from the surface of 
its primary, the energy increment due to propulsion is indicated 
in the display. Similarly, the decrement in orbital energy ef- 
fected by rocket braking is indicated. And, as a nonself-pro- 
pelled vehicle dips in and out of the atmosphere in a series of brak- 
ing ellipses, the energy indication drops at each perigee passage. 
Finally, when the altitude indication reaches 1 upon contact with 
the primary’s surface, energy is reduced to —1.0. At that alti- 
tude, the energy can exceed —1.0 only if the vehicle’s velocity is 
greater than 0.0, a state of affairs likely to be terminated abruptly. 


On the Viscous Flow Around the Leading Edge 
of a Flat Platet 


Bruno A. Boley and Morton B. Friedman 
Institute of Flight Structures, Columbia University, New York 
February 2, 1959 


gen FLOW in the neighborhood of the leading edge of a flat 
plate, even in the incompressible case and at zerc angle of at- 
tack, has been a problem of considerable theoretical interest ;' 2 
in addition, the importance of considering problems of this type 
as stepping stones toward further research has recently been 
pointed cut by Goldstein.’ The purpose of this note is to outline 
an analysis of the leading-edge flow regime, which though ap- 
proximate, appears to contain the essential features of the prob- 
lem. * 

Asa preliminary step, we shall consider the well-known Blasius 
solution from the viewpoint to be adopted in the leading-edge 
domain. The Prandtl equation will be satisfied if 


L h(x) 
f, f, [u(Ou/Ox) + v(Ou/Ov) — dydx = 
(1) 
for any arbitrary virtual variation 6u, where h(x) is the boundary- 
layer thickness, and L the length over which a solution is desired 
The appropriate boundary conditions, as well as the equation 
of continuity, are satisfied if the velocity components are derived 
from the stream function 
W(x, vy) = Uh/8[6(y/h)? — (y/h)4} (2) 
where LU’ is the free-stream velocity. 
Substitution of this function in Eq. (1) with 6u = (07y/(Oydh) X 
6h, and integration with respect to y, give 


L 


where primes indicate differentiation with respect to x. Eq. (3) is 
satisfied by the selutien of 


hh’ = 64v/7U (3a) 
which is 

h = k(vx/U)"? (3b) 
with k = 4.275 (3c) 


The corresponding displacement and momentum thicknesses 
are, respectively, 6* = 1.60V yx/U and 6 = 0.60V vx /U which 
compare favorably with the corresponding Blasius values? which 
are 6* = 1.73 vx, U and 0 = 0.66V vx/U. 

As an alternate procedure, we may assume a form for the func- 
tion A(x) in terms of some arbitrary parameters, and carry out the 
integration indicated in Eq. (3). If his chosen as that of Eq. (3b), 
with k an unknown parameter, then the same value for k, and 


+ This work was supported by the National Science Foundation under 
grant number NSF-G5948 
* The authors acknowledge the valuable suggestions of D. Dicker during 
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therefore the same solution as above, will, of course, result. The 
method is relatively insensitive to different chcices in the assumed 
function, as may be verified by performing the calculations on the 
basis of a different function h. For example, if we arbitrarily select 


h(x) = (4) 
with @ the required parameter, then the result is 
a = (da) 


and, therefore, depends here (as it did not in the previous case) on 
the length Z. The results have meaning only for x < L, where 
they are in good agreement, except near the leading edge, with 
the ones of Eqs. (3b), (c); this can be seen by the comparison of 
the boundary-layer thicknesses, obtained by the two methods, 
plotted in Fig. 1. Of course, for better accuracy for large L, a 
form for h more general than that cf Eq. (4)--i.e., with more 
parameters, needs to be chosen. The drag coefficients correspond- 
ing to the solutions of Eq. (3b), of Eq. (4a), and of Blasius are all 
of the same form and are, respectively, 1.4/VR1, 1.23/VR, and 
1.338/V where = UL/v. 

In the leading-edge regime, the Prandtl equations are not valid 
and give rise to singularities in the velocity field; the use of the 
complete Navier-Stokes equations can be expected” to eliminate 
this type of behavior.f We start by assuming that the boundary 
layer extends a distance a (as yet unknown) ahead of the leading 
edge along the axis of the plate (y = 0). Along y = 0, the velocity 
decreases continuously frem the free-stream value U at x = —a 
to zero at the leading edge x = 0. At the same time, the bound- 
ary layer grows from zere thickness at x = —a ina manner to be 
discussed. 

A simple flow field which satisfies these conditions is given by 
the stream function 


Uf ly (w + a)?/2ha?| + 


y — [Wx + a)*/a*]}; <0 

OS 
OSx<L 


t This problem was considered by Kuo.‘ His solution, obtained by a series 
approach combined with Lighthill's method for improving higher approxi- 
mations, shows the great complexity of the problem. Kuo’s results are in 
yeneral agreement with the ones described here. 


where h(x) is the boundary-layer thickness. The unknowns are 
the parameter a and the function A(x); asa first approximation } 
was taken to be 
h(x) = kV (x + a)(v/U); 

namely, that of Eq. (3b) shifted a distance a. The value of 
is taken as in Eq. (8c). The only parameter tc be determined 
therefore, is a. 

The terms of a stream function, the Navier-Stokes equations 
require that 


L(y) = — — 
— — 
(dy = 0 (6) 


This is clearly satisfied if 


Lp h(x) 
f L[pléy dydx = 0 (7) 
—-aJ0 


for any arbitrary virtual variation 6y. Substitution of Eq. (5) 
into expression (7) gives upon integration and simplification the 
following for R, = aU/v: 


+ + + + 1}*/2} 
where 


R, = LU/vand C, = 425/1872; C, = 2399/1055; 
Cs; = 153/400; C, = 63/80; C; = 19/48; Cs = 1/3 


The only admissible root for all chcices of R, is that plotted in 
Fig. 2, and is virtually independent of the choice of R, for R, > 
40. The boundary-layer height at x = 0 is given by hU/y = 
kV R,; away from the leading edge, the boundary layer over the 
plate very closely resembles that of the classical theory. As a 
consequence, the drag is almost unaffected by the introduction of 
the nose region. 
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Experimental Verification of Boundary-Layer 
orrections in Hypersonic Nozzles 


Donald L. Baradell 

Aeronautical Research Engineer, Langley Research Center, 
NASA, Langley Field, Va. 

February 2, 1959 


HE PROBLEM of accurately predicting boundary-layer growth 
for supersonic nozzle design assumes increased importance 
as the design Mach Number of the nozzles is pushed into the hy- 
personic regime. Several methods of calculating this growth have 
been advanced which vary in ease of computation and accuracy 
of results. Results of the calibrations performed on two hyper- 
sonic helium nozzles recently put into operation in the 11-in. 
Hypersonic Tunnel Section at the NASA Langley Research Cen- 
ter indicate that the method employed for boundary-layer calcu- 
lations in these nozzles is adequate. The flow in the test region 
of both nozzles has been found to be of good quality and the 
Mach Numbers obtained agree very satisfactorily with the de- 
sired design Mach Number. 
These nozzles were designed by the method of characteristics' 
to produce uniform parallel flow at Mach Numbers of 10 and 18, 
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the computations of the characteristic net being performed on an 
IBM 704 calculator. Both nozzles are axisymmetric with a 10.5- 
in. diameter at the center of the test region. The flow angles at 
the inflection points of tne nozzles were specified as 10° for the 
Mach-10 nozzle and 12° for the Mach-18 nozzle. The method of 
computing turbulent boundary-layer displacement effects devised 
by Persh and Lee? was used to obtain the boundary-layer displace- 
ment thickness in these nozzles after being adapted for high- 
speed machine computation. This method, which is based on a 
finite difference solution of the von Karman momentum equation, 
is applicable to either two-dimensional or axisymmetric nozzles 
and includes the effects of heat transfer. The two problems were 
set up so that after the ordinates of the flow field were obtained, 
these results could be used directly as input for the boundary-layer 
calculations. In order to obtain the desired test section dimen- 
sions, several iterations were necessary. 

As shown in Fig. I(a), the Mach-10 nozzle was designed to 
operate at a stagnation pressure (p,) of 100 psi. Because of pres- 
sure ratio difficulties the nozzle has not been calibrated at this 
pressure, but an extrapolation of the results obtained at higher 
pressures indicates that a Mach Number of 10 would be obtained 
if the nozzle were run at design pressure. The results of the cali- 
bration of the Mach-18 nozzle, as shown in Fig. 1(b), indicate 
that at the design pressure of 1,000 psi the average Mach Num- 
ber in the center of the test region is 17.8, which differs from de- 
sign Mach Number by only 1 per cent. 

In order to obtain a larger testing area, the length of the Mach- 
18 nozzle was limited by a method similar to that presented by 
Kenny and Yu.* The agreement between the design longitu- 
dinal Mach Number distribution for the Mach-18 nozzle and 
values obtained experimentally is shown in Fig. 2. The Mach 
Number gradient along the nozzle axis in the test region is only 
0.05 per in. in the Mach-18 nozzle and 0.01 per in. in the Mach-10 
nozzle. 

The pressure distribution in the test region of both nozzles was 
determined by a survey with an impact pressure rake. As shown 
in Fig. 3, the distribution in both nozzles is good, and a moderate 
change in stagnation pressure effects only a overall shift in the 
level of the distribution. 
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Figure 1. —Comparison of Mach numbers obteined by 
impact pressure surveys with design Mack numbers 
for two hypersonic nozzles. 
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Figure 2. — Longitudinal Mach number distribution 
in the Mach 18 helium nozzle. 
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Figure 3. -Results of & survey with ean impact pressure reke in 
the center of the test region in two hypersonic helium nozzles. 


Static pressures obtained along the nozzle wall during the cali- 
bration runs agree well with the design wall pressures and sub- 
stantiate the results obtained with the impact pressure rake. The 
values of Mach Number, J/, in the figures were obtained from the 
impact pressure ratios and are true values in all cases considered 
here, except those where the boundary layer influences the 
measured impact pressure. Accordingly, in Fig. 3(b), the Mach 
Number scale applies only to a region within about 2 in. of the 
nozzle axis of symmetry. 
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Remarks on ‘‘The Solution of the Laminar 
Boundary-Layer Equations”’ 


M. Z. v. Krzywoblocki 
Professor, University of Illinois, Urbana, III. 
February 13, 1959 


| era DEVELOPED a new iteration procedure for solving the 
so-called Prandtl boundary-layer equation. I would like to 
add a few remarks on this subject and on the subject of the 
boundary-layer equations in general. 

By means of a not one-to-one transformation of coordinates, 
Prandtl and Blasius passed from a partial differential equation to 
an ordinary one. All the procedures referring to a solution of the 
Prandtl-Blasius equation (Epstein’s as well) refer to an ordinary 
differential equation. Due to the existence of complex boundary 
conditions and the lack of a one-to-one transformation, it is not 
possible—using the available methods of the theory of functions 
to prove anything back in the domain of the partial differential 
equation. This can be done by means of some algebraic methods 
(Michal,? Morgan,’ ef al.) but only for an equation without any 
boundary conditions. Thus, an exact solution in the domain of 
an ordinary differential equation does not need to be necessarily 
an exact solution of the original partial differential equation. As 
a matter of fact, absolutely nothing can be said on that at the 
present time. There are some voices here and there, that the 
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future development of topology may give some answer, if any, 
on that question. 

The second remark refers to the Blasius technique. Blasius at- 
tempted to solve the ordinary differential equation of the bound- 
ary layer by matching together at an appropriate point a power 
series development, usually around the origin, with an asymptotic 
development around the point at infinity. But already in 1942 
H. Weyl* ® had shown that this technique may be erroneous (in a 
few cases he had shown that it 7s erroneous) since those two series 
may never meet together (a gap between them). Anyhow, Blas- 
ius proposed only a construction of a solution without providing 
any proofs of existence of a solution and possibly of uniqueness, 
whatsoever. 

It was H. Weyl who first attacked these more fundamental 
problems of the mathematical nature of the boundary layer 
and showed that, instead of using Blasius’ technique, it is better 
to use his own alternating series, which is very rapidly convergent. 
It seems to me that aerodynamists should be acquainted with 
Weyl’s alternating series since it is definitely superior to the 
Blasius technique. 

The third remark refers to the convergence of the iteration 
procedure of Epstein and to his remark on the subject of its con- 
vergence. Under suitably chosen conditions and bounds, it 
should be possible to prove—at least, a local-—convergence of the 
procedure in question. 

All the remarks given above are elaborated very thoroughly in 
many papers of mine on the mathematical aspects of the bound- 
ary-layer theory.6~"” 
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On a Simple Relation of Exact Airfoil Theory 


William T. Evans 

Aeronautical Research Scientist, Ames Research Center, NASA, 
Moffett Field, Calif. 

February 16, 1959 


XACT THEORY for incompressible potential flow yields the fol- 
lowing relation for the distribution of surface velocities about 
a two-dimensional airfoil 


U = k{sin (¢ + a) + sin (a + B)] (1) 
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where LU is the (transformed) velocity made dimensionless with 
respect to free stream, k is a function of airfoil geometry only, ¢ 
the angular coordinate on a circle conformally related to the air- 
foil, a the angle of attack, and 6 the negative of the angle of zero 
lift. By simple trigonometry, an interesting relation can be 
derived which has apparently been overlooked, and which bears 
on concepts arising from thin-airfoil theory. If the notation 
a* =a + B, o* = ¢ — Bis introduced, then by expanding and 
collecting terms, 
U = k[sin cos a* + (1 + cos ¢*) sin (2) 
so that at zero lift (a* = 0), 
Uo = k sin ¢* (3) 
and at a* = 7/2, 
U; = k(1 + cos ¢*) (4) 
(The velocity distribution U; corresponds to the maximum dimen- 
sionless circulation to which the Kutta condition can give rise.) 
Evidently, then, 
U = Uy cos a* + U; sin a* . (5a) 
and since lift coefficient is exactly proportional to sin a*, 
U = a* + (5b) 


where the final subscript 1 again refers to a* = 7/2. 

It is clear that, for small a*, the two terms of either Eq. (5a) 
or Eq. (5b) are somewhat analogous to the thin-airfoil concepts 
of ‘‘basie’”’ and ‘“‘additional’’ velocity distributions, respectively. 
In fact, for symmetrical airfoils, the factor U\/c;, is the same as the 
familiar ‘‘additional velocity increment’ Av,/V of reference 3. 
The correctness of Eq. (5b) for such airfoils has already been 
demonstrated by Loftin,* using the notation of reference 3, al 
though Av,/V was not interpreted as U;/c;,. For cambered air- 
foils, of course, these two quantities are not the same, since Av,/V 
is computed for uncambered thickness forms only, in accord with 
thin-airfoil concepts, and can no longer be used in an exact man- 
ner for cambered airfoils. Eqs. (5), however, are valid for all air- 
foils. 

In thin-airfoil theory, the basic velocity distribution (that to 
which the ‘‘additional” is added) is chosen in an essentially arbi- 
trary manner, usually at the “ideal’’ angle of attack® rather than 
at zero lift... A somewhat similar arbitrariness can be demon- 
strated in exact theory by generalizing Eq. (5a). While the 
generalized equation is easily derived from Eq. (1), it may be in- 
structive to use a different approach. 

Let the ‘‘chord” of an arbitrary airfoil be a completely arbi- 
trary reference line—that is, the flow at a = 0 is any arbitrarily 
chosen flow which satisfies the Kutta condition, with surface- 
velocity distribution U,. At @ = 7/2, the “chord” is normal to 
free stream, and the velocity distribution may be designated U’,. 
At arbitrary a, then, the free-stream velocity vector may be re- 
solved into components parallel and normal to the ‘‘chord,”’ 
and the distributions due to each component may be added. The 
result is 


U = U,cosa+ U, sina (6) 


For small a, the concepts of ‘basic’ and ‘additional’ velocity 
distributions again have an approximate validity for the two terms 
of this equation. Only if the ‘‘basic’’ c; is small, however, will the 
second term be approximately proportional to ‘‘additional”’ c;. 
Furthermore, the concept that the effect of airfoil shape (either 
thickness, camber, or both) is primarily embodied in the ‘“‘basic’”’ 
velocity distribution finds no interpretation in the equaticns 
above. 

Certain further concepts are suggested by the equations in their 
own right, and these may be of interest. The ratio ci/c:, in Eq. 
(5b) is the ratio of any given lift coefficient to the maximum to 
which the Kutta condition can give rise. The term ‘circulation 
ratio”? and symbol «x seem appropriate for this quantity. In 
terms of x, Eq. (5b) becomes 


U 2 + Vx (5c) 
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The concept of circulation ratio might be useful in applications 
of boundary-layer control and forced circulation, in which ratios 
approaching unity, and effective ratios exceeding unity, respec- 
tively, might be encountered. The use of x would lend emphasis 
to the theoretical nonlinearity and low slopes of lift curves for 
such flows, since x = sin a* (when the Kutta condition is satis- 
fied ) 

The use of the coordinate g* in Eqs. (2)—-(4) serves to emphasize 
the unique conformal correspondence among all airfoils which 
maps trailing edges into trailing edges, and which includes the 
circle considered as an airfoil with “trailing edge’ at ¢* = 7. A 
characteristic of this conformal correspondence may be illustrated 
by considering the division of Eq. (4) by Eq. (3): 

U,/Uo = cot(¢*/2) (7) 


This ratio is the same for all airfoils; furthermore, the ratio b<- 
tween the velocity distributions for any two absolute angles of 
attack a* is likewise independent of airfoil geometry, as can be 
seen from Eq. (2). This means that, in principle, the conformal 
correspondence between any airfoil and a circle, and, therefore, 
between any two airfoils—e.g., the same airfoil with different flap 
deflections—can be established approximately from measured 
pressure distributions for two different angles of attack. From 
the same data, computations of pressure distribution for other 
angles of attack could also be carried out. 

Finally, it may be of interest to consider the use of « and a*, 
Uy) and U; in computing flows in which both incidence and lift 
are specified—e.g., to agree with experiment—the Kutta condi- 
tion and results near the trailing edge being ignored. This is done 
by replacing sin(a + 8) in Eq. (1) by the required circulation, 
made dimensionless with respect to the circulation for coincident 
stagnation points on the circle;! this is simply «x: 


U = k{sin(g + a) + «] (8) 

From this, it is readily shown that 
U = (Uy cos a* + U; sin a*) — R(sin a*® — x) (9) 
which gives the result as the sum of a Kutta condition term and a 


correction term. 
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An Investigation of the Pressure Distribution 
on Conical Bodies in Hypersonic Flows* 


Victor Zakkay 

Aerodynamics Laboratory, Polytechnic Institute of Brooklyn, 
Freeport, N.Y. 

February 9, 1959 


LARGE AMOUNT of work on conical flow fields without axial 
symmetry at supersonic speed is presently available. How- 

ever, no apparent hypersonic approximation has yet been derived. 
* This research was supported by the USAF through the Air Force Office of 
Scientific Research, Air Research and Development Command, under 
Contract AF 49(638)-217, Project No. 9781, under the guidance of Dr. 
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Fic. 2. Pressure distributions on an ellipitical cone at a = 0°, 
10°, and 20° at M.. = 6 


In this note, experimental data on two elliptical cones at MV, = 6 
are presented and a hypersonic approach obtained from physical 
considerations is suggested. 

The method consists of determining at each point the local 
normal component of the free-stream velocity with respect to 
the surface of the body. Then, an equivalent circular cone at 
zero angle of attack is determined that has the same normal com- 
ponent of the free-stream velocity. The local pressure coefficient 
may then be obtained from the data available for circular cones 
at zero angle of attack and the same Mach Number. The re- 
sulting pressure distribution, as determined by this method, 
compares more favorably with the experimental results than the 
Newtonian-theory presented in reference 1. 

The pressure distributions on two elliptical cones having a 


* 16.7° equivalent circular cone, and ratios of the major to the 


minor axis of a/b equal to 1.39 and 1.87, are presented in Figs. 1 
and 2, for angles of attack of 0, 10°, and 20°. The experiments 
have been perforined at the Polytechnic Institute of Brooklyn 
Aerodynamics Laboratory. The ratio of the height of the cones 
to the square root of the base area for both models was H/( rab)'/* 
= 1.87, where H is the height of the cone. The stagnation con- 
ditions in the tunnel were 615 psia and approximately 1,700°R. 
for all tests. The Reynolds Number per in. based on free-stream 
conditions in the test was 0.31 & 10°. A summary of these re- 
sults has been presented in reference 2 and the detailed results 
of this investigation have been reported in reference 3. 

The results of the experiments are compared with the New- 
tonian theory and with the method presented here. 
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Comment on Nonlinear Pressure Method 
for Unsteady Aerodynamic Forces 


Harry L. Runyan and Homer G. Morgan 
Head, Dynamic Loads Analysis Branch, and Aeronautical Research 


Engineer, Respectively, 
Langley Research Center, NASA, Langley Field, Va. 


February 13, 1959 


N REFERENCE 1, several procedures for computing the air 
forces on an oscillating two-dimensional airfoil at high flight 
speeds were given. One method was termed the nonlinear pres- 
sure method. This method consisted of the insertion of the 
linear velocity potential in certain terms of the nonlinear pres- 
sure coefficient. Satisfactory agreement with other theories was 
shown by use of the method in flutter calculations. It was specu- 
lated and hoped that the method might be extended to the three- 
dimensional case and would provide a method for accounting for 
both airfoil thickness effects and finite aspect ratio. However, 
upon attempted application to the three-dimensional case, cer- 
tain inconsistencies were noted and were called to the attention 
of the authors by Garabed Zartarian, of Massachusetts Institute 
of Technology. It is now felt that the use of this method in the 
three-dimensional case, as indicated, is not warranted. 
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A Total-Temperature Probe for 
High-Temperature Boundary-Layer 
Measurements 


M. Sibulkin 
Convair Scientific Research Laboratory, San Diego, Calif. 
February 20, 1959 


N ORDER TO OBTAIN an exact measurement of the total tem- 
perature of a moving fluid, the fluid would have to be brought 
to rest adiabatically and its temperature then measured by a 
thermosensing element which did not affect the temperature of the 
fluid. It has not been possible to devise an instrument to do this, 
and, in practice, total-temperature probes suffer from errors due 
to nonzero fluid velocities past the sensing element, radiation of 
heat from the probe to its surroundings (important primarily at 
elevated fluid temperatures), and conduction of heat along the 
sensing element.'! For boundary-layer surveys, the requirement 
that the probe dimensions be small compared to the boundary- 
layer thickness aggravates the problem, and it is necessary to 
calibrate the probe against known total temperatures. 

In order to reduce these calibration corrections, an attempt has 
been made to design a total-temperature probe for use in super- 
sonic flow (Fig. 1) which minimizes the errors mentioned above. 
The probe is essentially a two-dimensional supersonic diffuser 
having a variable exit area (controlled by translating the exit 
plug). After the velocity of the air has been reduced to a low value 
(M< 1) by the diffuser, its temperature is sensed by the thermo- 
couple glued to the mica crosspiece. The thermocouple has been 
laid in the manner shown, and mica has been used to reduce the 
conduction of heat to probe frame. The two-dimensional shape 
was chosen (for use in two-dimensional boundary layers) to mini- 
mize the probe dimension in the direction perpendicular to the 
wall while still allowing use of the relatively long thermocouple 
shown.* An earlier probe,? constructed for use in the heated air- 


* The two-dimensional geometry has been utilized independently by 
Bradfield‘ for low-temperature boundary-layer measurements. 
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stream of the Naval Ordnance Laboratory hypersonic wind tun- 
nel, was goldplated to minimize radiation errors. 


CALIBRATION MEASUREMENTS 


Measurements of the probe temperature 7, were made in the 
uniform flow section of the Jet Propulsion Laboratory 20-in. 
supersonic wind tunnel over a Mach Number range M = 2.3 to 
5.0 at the available total temperature 7, = 100 + 20°F. It was 
assumed that 7; was equal to the temperature measured in the 
low-speed flow at the tunnel supply section. The results are 
shown in Fig. 2 where the temperature ratio 7,/7; is correlated 
with the density inside the probe p, for reasons discussed later, 
For a given Mach Number, p, was varied by changing the tunnel 
supply pressure. The air velocity past the sensing element was 
varied by changing the exit height ,—i.e., the distance between 
the exit plug and the side plates; the exit area A, = 2 X 0.5 X h, 
in.2 Increasing the internal air velocity (1) increases the rate of 
heat transfer to the sensing element which reduces conduction 
and radiation errors,! and (2) decreases the internal air static tem- 
perature. In general, 7/7, will have a maximum at some inter- 
mediate air velocity. The data given in Fig. 2 are for the two 
values of h, nearest this maximum. 

When the external flow is supersonic, the flow inside the probe is 
choked—i.e, the Mach Number at the probe exit 1/7, = 1. There- 
fore, conditions inside the probe are only indirectly related to the 
external stream Mach Number and Reynolds Number. However, 
for a given total temperature, it is the Reynolds Number based 
upon conditions inside the probe Re, which determines the rate of 
heat transfer to the sensing element and, consequently, 7). De- 
fining 


Rey = (ppttyl)/mp = 


where a is the speed of sound, yu the viscosity, / a characteristic 
length, and the subscript refers to conditions in the vicinity of 
the sensing element; the following estimate of Re, may be made. 
Since M, = 1, 


M, = f(A,/Ae) 
Also 
ap = (yRT,)!/2 


and assuming a power law variation of viscosity 


Fic. 1. Photograph, tip details, and exploded view of the total- 
temperature probe. 
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Fic. 2. Effect of air density inside the probe p, on measured 


total temperature 7). 


Mp = 
Substituting these relations in Eq. (1) and knowing that 7, ~ 7, 
gives 
Rep = ~ “py (2) 


For a fixed level of 7;, Eq. (2) gives Re, ~ pp» which suggests the 
use of p, to correlate the data in Fig. 2. For significantly higher 
levels of T;, lower values of 7/7; would be expected. 

Before using such a probe for boundary-layer measurements, its 
interference effects (as discussed, for example, in reference 3) 
should be considered; such effects are known to be larger for 
laminar than for turbulent boundary lavers. 
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Eddy Diffusivities in Forced and Free 
Convection Boundary Layers 


D. E. Bourne 

Assistant Lecturer in Applied Mathematics, University of Sheffield, 
Sheftield 10, England 

February 13, 1959 


HE PURPOSE of this note is that of discussing the widely used 
g penencoren that in incompressible turbulent boundary layers 
over heated surfaces the eddy momentum diffusivity, e, is identi- 
cal in value with the eddy heat diffusivity, e7. This assumption 
appears to have arisen in the following way. 

In 1874, Osborne Reynolds! suggested that the rate at which 
heat is transferred to a fluid flowing over a heated surface is 
directly proportional to the friction between the fluid and the 
surface; this assumption is ‘Reynolds’ Analogy.’”’ In the case of 


heat transferred by forced convection through an incompressible 
turbulent boundary layer over a plane surface (two-dimensional 
problem), the analogy is valid? if (a) the mean downstream pres- 
sure gradient is zero, (b) the Prandtl Number, ao, is unity, (c) 
the boundary conditions for the mean temperature, 7, and the 
downstream component of mean velocity, U, are similar and 


(d) « = ev. (Here the concept of eddy diffusivities has, of course, 
been assumed, and the momentum and temperature equations 
are used in the forms given by the momentum-transfer theory ; 
for practical purposes, both these assumptions are now generally 
accepted. ) 

Since Reynolds’ Analogy has been used with satisfactory re- 
sults in a number of cases, it is reasonable to inquire whether the 
condition ¢ = ey holds generally. 

In the case of forced convection through an air (¢ = 0.7) 
boundary layer on a flat plate, Davies* supposed that ey = Xe, 
where A is some constant. It was then shown that agreement be- 
tween theory and experiment for the temperature distribution in 
the boundary layer was achieved by taking A = 1. This strongly 
suggests that in flows under similar conditions we may suppose 
= 

However, there are certain cases in which the assumption can- 
not hold. Consider, fer example, the flow over a plane surface. 
Taking axes Ox, along the surface in the direction parallel to the 
mainstream, and Oy, normal to the surface, we define 

e = —u'v'/(OU/Oy) ex = —v'T'/(OT/dy) 
(u’, v') are the eddy velocity fluctuations about the mean com- 
ponents (U, V), and 7” is the eddy fluctuation of the temperature 
about the mean value, 7. 

Now, in the flow over a heated disc rotating in a fluid at rest at 
large distances from it, if we suppose x to be the radial distance 


from the center, ’ = Oat y = Oandat y = ~. Thus, OU’/dy 
vanishes at some point, and since 07'/Oy vanishes only at y = &, 
the condition « = ey will break down at (and near to) this 
point. 


It has also been assumed by some writers—e.g., Bayley*-that 
in the free-convection boundary layer on a vertical heated 
flat plate in air we may equate the eddy momentum and heat 
diffusivities. However, with coordinates as before, so that x 
measures distance along the plate, 0LU’/Oy vanishes at some points 
in the boundary layer, and tke condition « = ey again breaks 
down. 

This free-convection problem has also been investigated by the 
present writer, using Eckert and Jackson’s® analysis of some ex- 
perimental data provided by Griffiths and Davis. Substituting 
their empirical formulas for the temperature and velocity dis- 
tributions, the mean momentum and temperature equations were 
integrated to determine ¢ and ey. The result shows that the 
value of en/\s varies from zero up to a maximum of about 5.5 
in the inner 50 per cent of the boundary tayer. The empirical 
formulas do not represent the experimental data sufficiently ac- 
curately in the outer 50 per cent of the layer to continue the 
analysis in this region with any confidence. It should also be 
emphasized that there is considerable scatter of the experimental 
results, so that the conclusions arising from them must be re- 
garded as tentative. (The data of the experimental work of 
Griffiths and Davis, performed in 1922, appears to be all that are 
available; further experimental werk in this field might be 
profitable. ) 

We conclude, therefore, that the assumption that the eddy heat 
and momentum diffusivities are identical in value should be used 
with considerable care. It is suggested, tentatively, that in 
problems of forced convection over plane surfaces the assumption 
may be used with some confidence. It is also not unreasonable to 
suppose that this result may extend to problems of forced convec- 
tion over transversely curved surfaces—e.g., forced convection 
through an axisymmetrical turbulent boundary layer on a circular 
cylinder. However, in cases in which the boundary conditions 
for the temperature and the downstream component of mean 
velocity are not similar, the assumption does not appear to be 
permissible. 
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Deduction of Mass-Mean Properties From a 
Knowledge of Pressure and Mean Density 


H. G. Elrod, Jr. 

Department of Mechanical Engineering, Columbia University, 
New York 

February 20, 1959 


be FOLLOWING PROBLEM has iiterest because of the possibility 
of using pressure volume measuremeiuts to determine high- 
temperature gas properties. Given a container of nonideal gas 
at uniform pressure, p, and nonuniform temperature, under 
what conditions can a knowledge of the pressure and mean 
density be used to deduce the correct mass-mean energy? 

Let us assume that the internal energy of the gas is expressed in 
the form 


u = f(p, p) (1) 


then an energy u* can be calculated from Eq. (1) using the pres- 
sure, p, and the mean density defined by: 


fdV/V (2) 
where V is the container volume. 
Thus: 
u* = f(p, p) (3) 
On the other hand, 
i= f'updV/pV (4) 
Equating u* and #, we obtain from Eqs. (3) and (4): 
= p)pdV (5) 
or: 
S — Mp, = 0 (6) 
Suppose, for the present, that 
p=ptp’ (7) 
where p’ is a small density deviation. 
Then: 


= f(p, 6) + (Of/Op)pp’ + + . . (8) 
Eq. (6) becomes: 
S + + dV = 0 (9) 
Now S =0 (10) 
Therefore, Eq. (9) reduces to: 
+ + O(p’%).. .JdV = 0 (11) 


Since p”? is positive and the derivatives in Eq. (10) are uniform 
throughout the volume, being evaluated at the state (/, ), it is 
necessary for 


(Of/Op), + (p/2)(07f/Op?), = O (12) 
By integration: 
u = f(p, p) = A(p) + B(p)/p (13) 


where A(p) and B(p) may be arbitrary functions. 
The equation of state, given by Eq. (13), is a necessary condi- 
tion for u* to equal #. Direct substitution of Eq. (13) into Eq. 
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(6) shows the condition also to be sufficient, whether or not p’ is 
small. Eq. (13) is a suitable caloric equation of state for a perfect 
gas having constant specific heats and for a number of noiideal 
gases and vapors. For all these, the mass-mean internal energy 
can be found from a knowledge only of the pressure and mass- 
mean density, without information on the temperature distribu- 
tion. In conclusion, it should be noted that the original question 
can be phrased in terms of extensive variables other than the 
internal energy. 


On the Normal Shock Wave Attached to the 
Curved Surface 


Takeo Sakurai * 

Department of Aeronautical Engineering, Faculty of Engineering, 
Kyoto University, Kyoto, Japan 

January 23, 1959 


yeaa ZIEREP! published a paper on the normal shock 
wave attached to a curved surface in two-dimensional flow 
of compressible inviscid fluid [Fig. l(a)]. He used orthogonal 
curvilinear coordinates of which the coordinate lines coincide 
with the shock wave and the body surface, respectively. He dis- 
cussed the shape of the shock wave near the surface by means of 
a limiting process, and found that the curvature of the shock wave 
must become infinite on approaching the surface. He also dis- 
cussed the fact that the configuration with the normal shock 
wave attached to the convex surface can exist in the Mach 
Number region 17 = 1 ~ 1.662. On the other hand, the same 
problem had been studied a few years before by Lin and Rubinoy,? 
using a somewhat different method, and analogous results had 
been obtained. They had concluded, however, that the con- 
figuration with normal shock wave attached to the convex surface 
cannot exist in the above Mach Number region, 7 = 1 ~ 1.662, 
in contradiction with Zierep’s results. Moreover, Zierep’s results 
on the entropy derivative and on several other derivatives behind 
the shock wave seem to be difficult to understand from the view 
point of common sense gasdynamics. The problem being of 
fundamental importance, it may be worth while to point out the 
cause of discrepancy between the two theories. 

In terms of the curvilinear coordinates, which were used by 
Zierep |Fig. 1(a)], the fundamental equations and the boundary 
conditions may be written as follows: 


+ — + + + 
h.WAp,/p) =0 (1) 


hWiWie + + + + 


(p/C,)St] = 0 (2) 
Win + + Sy + + 
(p/C,)Sn] = 0 (3) 


* The author wants to express his cordial thanks to Prof. Ko Tamada and 
Prof. Isao Imai for their guidance and encouragement. He is also indebted 
to Prof. Hidenori Hasimoto for his invaluable advice and discussion. 


Fic. l(a). Fic. 1(b). 
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hy Wee + h, + he W By = —hi( T/W2)Sé | 
h(T/Wi)S_ 


p = RpT, p = (5) 
=F on&£=0 (6) 

+ (2/2) = (1942/2) 
W. =0...onyn = 0 (7) 


From the above equations, we can determine the flow conditions 
behind the shock wave—i.e., 1, Ws, ete., provided that the shape 
of the shock wave together with the flow conditions in front of the 
shock wave are given. Therefore, by differentiating these along 
the shock wave (£ = 0), we may obtain the y-derivatives of Wi, 
W», etc., immediately behind the shock. The £-derivatives can 
be found from the y-derivatives just obtained and the basic equa- 
tions. Then the question is whether the curvature of the stream- 
line behind the shock wave may coincide with that of the surface 
itself in the limit of approaching the surface, or not. (The curva- 
ture of the surface is assumed to be continuous in passing the 
shock wave.) 

Now, the difference between the curvature of the streamline 
and the coordinate line 7 = const is given as follows [Fig. 1 (b)]: 
6 = = — Weg) /( + We?) (8) 
This 6 should vanish in the limit 7 — 0 in order that the fluid may 
follow the body surface. Thus, we obtain the condition: 

(Withee — = 0 (9) 

0 
from which the shape of the shock wave near its root may be de- 
termined. The condition (9) coincides with that of Zierep only if 
W.Wi¢ tends to zero in the said limit. As was shown by Zierep 
himself, however, W, Wis, unfortunately, tends to a finite value. 
Hence, we cannot admit Zierep’s condition Wz — 0; instead, we 
must adopt Eq. (9) as it is. 

If we correct Zierep’s analysis, using Eq. (9) in place of his 
erroneous condition, we arrive at the following results: 

hoW = —G&h, (10) 


(1/C,)heSy = (11) 


where & and 7 are the functions of 1/ only, which are shown in 
Fig. 2. From the above results, we can proceed in the same way 
as Zierep and find that the configuration with normal shock wave 
attached to the convex surface cannot exist in the Mach Number 
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region J = 1 ~ 1.662, in complete agreement with the results of 
Lin and Rubinov. 
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A Mechanism of Resonance Tubes* 


J. Wilson and E. L. Resler, Jr. 
Corneli University, Ithaca, N.Y. 
February 20, 1959 


SYMBOLS 
‘ = speed of sound 
Cp = specific heat at constant pressure 
f = frequency 
F(t) = velocity as a function of — at ¢ = 0 
G(t) = sonic velocity as a function of — att = 0 
l = length of resonance tube 
p = pressure 
R = gas constant 
S = entropy 
7 = temperature 
t = time 
u = velocity 
¥ = ratio of specific heats 
é = distance along tube 
Subscript 
0 = average 


INTRODUCTION 


Tre RECENT PAPERS": 2 suggest that the mechanism of heating 

inside a resonance tube may be explained by the passage of 
shock waves along the tube. The photographs of reference 2 
appear to support this claim. Similar pictures have also been 
obtained at Cornell. It is the purpose of this paper to show 
that a criterion of wave steepening puts an upper limit to the 
temperature attainable in the tube. 

The photographs of the waves in a resonance tube obtained by 
Hall and Berry seem to confirm the assumption of steepening 
pressure waves rather than shock waves, as the disturbance grows 
narrower on its passage down the tube, as would be expected for 
a steepening pressure wave. It is not clear how this could be 
explained if the disturbance were a shock. 


WAVE STEEPENING 
Consider a closed tube with the pressure at the open end 
varying periodically: 
p = po + Ap(t) 
An increase of entropy of the air in the tube is expressed by 
dS = C,(dT/T) — R(dp/p) 


and for any number of complete cycles S (dp /p) = 0, so 
dT/T = dS/C,, 
putting dS = ASfdt where AS is the entropy increase per cycle. 
dT/dt = (TAS/C,)f 
Thus, the temperature will rise until the entropy increase per 
cycle becomes zero. The entropy increase is caused by the pres- 
sure waves inside the hole having steepened to become shock 
waves, and so the temperature rises until the waves can no longer 


* This work was carried out at the Graduate School of Aeronautical 
Engineering, under the sponsorship of the Mechanics Branch, Office of 
Naval Research. Reproduction in whole or in part is permitted for any 
purpose of the United States Government. 
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steepen sufficiently to become shock waves within the length of 
the tube. 

From Eq. (41.04) of reference 4, steepening occurs at time, é, 
such that 

1 + [F(é) + G(é)]t = 0 
For steepening at the exit, 
t = 2l/(u + c) = 21/|F(é) + G(é)] 

Setting = 0 in Eqs. (40.10) and approximating 

F(t) + = + = col + + 1)/2y] [Ap(t)/pol} 
Thus, 
+ GE) = + 1)/2y]- Col (d/dt) [Ap(t)/pol} -(dt/dé) 
+ 1)/2y]-Cof(d/dt) [Ap(t)/pol} X 

{1/[F(E) + 


and the steepening condition becomes 
[2ley(d/dt) [Ap(t)/po] | 
1)/2 = 0 
With co? = yRT>, this gives 
To = (1/yR) + [Ap(t)/po] + 
}1 + + 1)/27] [Ap(t)/pol } 2) max? = 
(1/yR) 1)/y] {(d/dt) [Ap(t)/pol} maxe 


This formula represents the maximum attainable temperature 
within the tube. If heat transfer to the sides of the tube occurs, 
then, in the final state, a finite amount of energy must be given 
to the air in the tube each cycle to replace energy lost by heat 
transfer. In this case, the waves will steepen to form shocks 
within the tube, and a lower final temperature will result. 

The formula indicates that there is a minimum rate of pressure 
increase which will cause heating. This is the value for which 
the waves steepen to a shock at the exit with room temperature 
in the tube. For lower rates of pressure increase, the waves 
cannot steepen to a shock within the hole and, thus, cannot cause 
heating. 


THE RESONANCE TUBE CYCLE 


The value of the rate of increase of pressure which is to be 
inserted into the temperature formula will be determined by the 
coupling between the internal and external flows. Fig. 1 shows 
a possible cycle for a resonance tube held in a sonic or supersonic 
jet. It is known that resonance takes place when the shock wave 
in front of the tube is in a region of instability. The steepened 
pressure wave emerging from the tube sets this shock in motion, 
causing it to move across the region of instability. This motion 
of the shock wave causes a decrease of pressure at the entrance to 
the tube, and so, an expansion wave travels down the tube. On 
returning, this expansion wave induces the shock to return 
through the region of instability, causing a compression wave to 
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travel down the tube. This compression wave on leaving the 
tube starts the cycle again. 

Thus, the maximum ratio of pressure increase [dAp(t)/dt] ma; 
depends on the motion of the shock wave under the influence of 
the expansion wave. A satisfactory calculation of this quantity 
has not been made. 


CONCLUDING REMARKS 


The above theory predicts a maximum temperature in a reso- 
nance tube and presents a possible coupling mechanism between 
the internal and external flows when the tube is held in a super- 
sonic jet. However, it does not explain the change of tempera- 
ture along the tube found by Sibulkin and Vrebalovich. Also, 
disturbances have been treated as small, which is not true for a 
resonance tube, and the effects of reflection at the end of the tube 
have been ignored. 

These factors and also the influence of viscosity should be taken 
into account in formulating a more complete theory. 
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On the Torsion of a Cylindrical Beam Having 
a Trapezoidal Cross Section 
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SUMMARY 


The Galerkin method is applied in the solution of a pure tor- 
sion problem. The stress field, displacement field, and torsional 
rigidity are determined for a cylindrical beam having a trape- 
zoidal cross section. A triangular cross section is considered as a 
special case. 

SOLUTION 


co SOLUTION of the problem will be approached through 
the Prandtl stress function using the Galerkin method. Then, 
for a beam as shown in Fig. 1, the stress field which satisfies the 
equation of equilibrium 


+ =0 (1) 
is given by: 
0 0 Gawy(x, y) 
y= 0 0 —Gavwr(x, y) (2) 
Gayy(x, ¥) —Gawz(x, 0 
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where G is the shear modulus, a is the twist per unit length, and Eq. (10). Then, Eqs. (7), (9), and (10) lead to: 
v(x, y) is the Prandtl stress function, the subscripts x, y denote e 
the corresponding derivatives, with respect to x and y. S Sir » Amnn(X, ») he af n(x, ydS = 0 (11) 
The displacement field is related to the W as follows: § m=1 
u = —ayz; v= axz; w= a®(x, y) (3) However, since 7»(x, y) are coordinate functions which have to 
aes satisfy the boundary conditions, Eq. (8), then Eq. (9) can be 
a written as: 
P(x, y) 
y) [Ov(x, y)]/dy} dx— { [dy(x, y)]/Ox} dy) (4) Wx, ¥) = nox, y)Pm(x, (12) 
and where n(T) = 0 (13) 
v(x, y) = W(x, y) + (1/2) (x? + 
Then, the torsional rigidity D is given by P(x, y) = Ai + Axx + Asyy +... + Anv?yt (14) 
Evidently 
D= 2cf f W(x, y)dS (6) 
m(x, ¥) = no(x, y) 
where S is the domain of integration. If W(x, y) is known, the 
solution of the problem is given by Eqs. (1)-(6). 
The Prandtl stress function is determined by finding the solu- 
tion for Poisson's differential equation nlx, y) = x? ytno(x, y) 
L[¥(x, »)] = V(x, 9) +2 = 0 (7) The Galerkin method, as indicated in the above section, will be 


applied to the case of a trapezoidal cross section. In this appli- 


in the region S, and subjected to the boundary conditions 
cation, because of limitation of space, we shall limit the selec- 


ar) = 0 (8) tion of W(x, y) to the first approximation, which is quite satis- 
on I’, the boundary of S. Hence, the solution of the problem factory for engineering purposes. However, no difficulties arise 
is determined as soon as V(x, y) is determined. A most conven- when taking two or more terms into account. Choosing » = 1, 
ient method is that of Galerkin" ?, using we have from Eq. (10) that 


n 
W(x, = Amnn(x, ¥) 9) 
(x, y) p> (m n (9) no(X, vas [T2no(x, ¥) lno(x, t (16) 
Ss Ss 


the Galerkin form of the stated problem becomes 5 ; : 
Note that for a trapezoidal section, Fig. 1, we have: 


ff L[W(x, y)]n(x, ydS = 0O(r = 1,2,3,...,) (10) n(x, = (cx — y)(y — dx) (a — x) (b — x) (17) 
S Then, 


where the n(x, y) are coordinate functions, which satisfy Eq. (8), 
and A,,’s are arbitrary constants which are to be determined from ~2 Sf no(x, y)dS = (c — d)3/180 [26° + a(3at> — 365 — 2a5)| 
s 


(18) 


ff IV 2no(x, ¥) |no(x, = 1/2520 — a) [952(c8 — + 475(c® d3) + 895(c3d? — + 15(c2d — cd*) + 
s 


4355(c'd — cd*)| + (ab? — a%b) [28(c8 — d®) + 776(c? — d3) — 845(c8d? — — 888(c%d — cd?) + 
4964(c'd — cd*)] + (a*b® — a®b?) [—266(c® — d5) — 14(c? — — 98(c3d? — c¥d3) + 42(c%d — cd*) — 196(c4d —cd*)]} (19) 


Eqs. (18) and (19) yield A; and, hence, %;. The value of the gi ee ve a 


torsi igidity is at termined by ing Eqs. (18) F ‘ ? ee 
( By setting a = 0 in the above equations, the solution is reduced 
and (6). Therefore, the stress field is given by : 3 : 

to the case of torsion of a beam of triangular cross section. 
= GaA;[On(x, y)/Oy] = GaA,|—2aby + Therefore: 


2(a + b)xy + ab(c + d)x + (c + d)x® — 
(a + b) (c + d)x? — 2x*y] 


—2 ff nox, y)dS | a=0 = (¢ — (22) 


24abcdxy — 18(c + d)x*y? + 8xy* — 36cd(a + 5) X 
+ 48cdx*y] + B 


ab(c + d)y + 3(c + d)xty — 475(c® — d®) + 895(c8d? — + 15(c8d — cd*) + 
2(a + + d)xy — 2xy*] 4355(ctd — cd*)] (23) 
The evaluation of the displacement field requires a more lengthy Then, for this case: 
calculation. Using Eqs. (5) and (9), Eq. (4) leads to: es = GaAy'[2bxy + (c + dx? — W(c + d)x® — 2x%y) (24) 
P(x, y) = A,/12[6ab(c + d)x? — 4a +b) (ec +d)X | + + + 
3(c + d)xty — 2b(c + d)xy — 2xy*] (25) 
| 


®'(x, y) = A,’/12[—4b(c + d)x? + 3(c + d)x* — 

4by? — 18(c + d)x*y? + — 36bcedx*y + 
where B is an arbitrary constant and corresponds to the rigid 48cdx*y] + B’ (26) 
body displacement in the Z-direction. The result from Eq. 
(16) when applied to Eq. (21) yields the displacement field to 
within rigid body motion as stated by the first boundary value 
problem in static elasticity. 


Note that the prime denotes corresponding values of stresses and 
displacement functions for the beams of triangular cross section. 
Other special cases can be handled by proper manipulation of the 


| 
tor- 
| 
ipe- 
AS a 
igh 
en, 
the 
(1) 
(2) 
JK 
= 
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four parameters a, 6, c, and d. For instance, when c or d = 0, 
and a = 0, we have the case of a right triangle. 
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= NOTE DEALS with an exact method for the calculation of 
local heat-transfer coefficients in a laminar boundary layer 
with constant properties, for arbitrary pressure and surface tem- 
perature distributions. The method depends on finding a 
formally exact power series representation for the temperature 
in terms of universal functions. Such a representation is possible 
if the exact solutions of the boundary-layer equations given by 
Gortler! and author? are used in solving the energy equation. 
As an illustration, the author’s solution will be used. 

The energy equation for incompressible steady flow over a two- 


dimensional body can be written as: 
pC, UT,. + VT, y) = UdP/dx + kT, y, (1) 
While the boundary conditions are given by, 
aty=0, T=T((x); asy>o, (2) 
The transformation given in reference 2 can be written as: 


= exp[—(1/r) fUa(x)dx], = yUol(x)/v, 
2) = u(x, y)/vy 
Letting 
T = T(x) + [To(x) — To(x)] OE, 7) (4) 
and using the above transformation, Eq. (1) reduces to 
with 
a(t) = —éd[In(7T) — o= pp, (6) 
where 7; is some reference temperature. The relation 
pC,dTo/dx = dP/dx (7) 
has been used in deriving Eq. (5). The boundary conditions for 
the function 6 are: 
0) = 1, O& ») as n> @ (8) 


The solution of Eq. (5) will be assumed as 


HE, = on )E" (9) 
0 


in the regime of convergence of 
0 
Substituting Eqs. (9) and (10) into Eq. (5), one finds 
6” — aogo’ = 0 
and 
On” + (m — a) = Qn-1, (n > 1) 


with 


JULY, 1959 


n—1 ™ 


Qn-1 = an—k + > am—k % + 
k=0 m=0 k=0 
n—1 


k=0 k=0 


The functions 6, satisfy the following boundary conditions: 
= 1,0(0) = 0 
and 
0,10) = = 0 for all | (14) 
For n = 1, the functions 6, can be expressed as linear combina- 
tions of universal functions. This is accomplished by letting 
= + Bihy 
= + + + Boh2 + 
= on + + asgs + + + 
+ + + aBirina + (15) 
where the §’s are the coefficients appearing in the expansion of 
A(é).2 If the differential operator L,,(z) is defined as 
L,(z) = + (n — ao)go'Z (16) 
the following differential equations for the functions g, h, r are 
obtained: 


n= Li(g1) = 
Li(hy) + 
= 2, 
L2(g2) (17) 
Lo( hy) + (ao — + 260'fir + fily’ 
Lo(he) + 26'fo 
Lo(r1,1) + + (ay —1 fi "gy + figs’ 


The boundary conditions for all the functions g, h, and r are 
g(0) = g~) = h(0) = = = = 0 (18) 
The class of functions 7)(x) T(x) considered in this 


analysis can be obtained from Eqs. (6) and (10) and the relation- 
ship connecting and x. It follows from Eqs. (6) and (10) that 


E~% ane" (19) 
0 


To 


but? 


hence, 


0 
where a,, bn, and C,, are constants. 
The heat-transfer rate at the wall is given by 


= —KT, = MTo — Ta) (22) 


where / is the local heat-transfer coefficient. Hence, using 
Eqs. (4), (9), and (22), one obtains 


h = —(koUa(x)/v) >, (23) 
0 


A procedure similar to the one given above can be applied to 
Gértler’s solution. Both solutions show that for 7) — To = f(x), 
the resulting universal functions depend on the two parameters, 
Bo, and ao, which characterize the free-stream velocity and the 
temperature difference, respectively. 


REFERENCES 


1 Gortler, H., A New Series for the Calculation of Laminar Boundary Layer 
Flows, J. Math. & Mech., Vol. 6, No. i, pp. 1-66, January, 1957. 

? Hassan, H. A., A New Solution io the Laminar Boundary-Layer Equa- 
tions, Readers’ Forum, Journal of the Aero/Space Sciences, Vol. 26, No. 3, 
pp. 189-190, March, 1959. 


fi 
— 
: 
2 
: 
: 
x! Bo b,x / Bo (20) 
0 
R E 
4 
be 
(10) 
(12) 


REQUIREMENTS for CONTRIBUTIONS 
to the publications of the 
INSTITUTE of the AERONAUTICAL SCIENCES 


The Institute of the Aeronautical Sciences invites both members and nonmembers from any 
country to submit papers for publication in the JouRNAL OF THE AgRO/Space ScrENCES and 
AERO/SPACE ENGINEERING. The Institute, following the practice of other societies, does not pay for 


contributions. 


The following directions for the preparation of papers, if followed by authors, will save corre- 
spondence, avoid the return of papers for changes, minimize the work of preparation for the printer, 
and save the expense due to the charges made for ‘‘author’s corrections.” 


Manuscripts: The original typewritten copy of the paper is 
desired. To expedite review, an additional copy of both the manu- 
script and figures should be submitted. The manuscript must 
be double or triple spaced on one side of white paper sheets, 
consecutively numbered. There should be wide margins to allow 
for the marking of directions to the printer. Correcting, chang- 
ing, or adding to papers after they are in type is costly. It is, 
therefore, imperative that papers submitted be in final form. 
Typographical errors may be corrected on proofs, but if authors 
wish to add material, they may do so at their own expense. In 
mailing, drawings may be rolled, but manuscripts should be sent 
flat. Send by first-class mail (register if you wish for your own 
protection) to the Editorial Office, Institute of the Aeronautical 
Sciences, 2 East 64th Street, New York 21, N.Y. All manuscripts 
will be examined by the Editorial Committee and by the Editor. 
Authors will be advised as promptly as possible whether the paper 
is acceptable for publication. 

Trries: The title of the paper should be brief. The name 
and initials of the author should be written as he prefers. The 
use of the full name of an author is advocated because of the fre- 
quent duplication of initials and surnames which sometimes makes 
it difficult to establish the identity of the author. This is par- 
ticularly important for large annual indexing and abstracting 
services. All titles and degrees or honors are omitted. The name 
of the organization with which the author is associated should 
be placed after his name on a separate line. The date on which 
the paper is received will be inserted by the Editor. The author’s 
title or position should be indicated in a footnote. 


SUMMARIES OR ApsTRACTs: An abstract to be printed at the 
beginning should accompany each article. It should be adequate 
asan index and asa summary. It should contain a statement of 
major conclusions reached, since summaries in many cases con- 
stitute the only source of information used in compiling scien- 
tific reference indexes. Abstracts printed in other journals, espe- 
cially foreign, in most cases consist of summaries from printed 
papers. The summary should explain as adequately as possible 
the major conclusions to a nonspecialist in the subject and should 
contain from 100 to 300 words, depending on the length of the 
paper. 

SuBHEADINGS: Subheadings should be inserted by the author 
at frequent intervals. The work of editorial preparation will be 
simplified by the author’s provision of many subheadings. 

MatTTER USUALLY DELETED: Photographs or illustrations of 
little technical interest and not showing advances in general 
practice. Too detailed tabular matter (general results of such 
tables may be included in the text). Lengthy descriptions of 
materials or processes or of preliminary experiments or theories 
that preceded final results; salient features only are of interest. 

REFERENCES: References should be numbered consecutively 
and grouped together at the end of the manuscript, with only 
the corresponding number being mentioned in the text. The 
arrangement should be as follows. For books: 'Durand, W. F., 
Aerodynamic Theory, 1st Ed., Vol. 1, p. 23; Julius Springer, 
Berlin, 1934. For magazines: ‘England, C. R., Crawford, 
A. B., and Mumford, W. W., Some Results of a Study of Ulitra- 
Short-Wave Transmission Phenomenon, Proc. IRE, Vol. 20, 
No. 12, pp. 481-482, March, 1933. Please give author, title, 
edition, volume, number, page, publisher, and place and date of 
publication as indicated. Omission of one required fact causes 
much extra editorial work and possible inaccuracies. 


ILLUSTRATIONS: Illustrations should accompany manuscripts, 
and each should always be referred to in the text by number. 
Drawings or graphs should not be larger than 12 X 16 inches, 
and must be made with jet black India ink on white paper or 
tracing cloth, the latter being preferred. Do not use typewriter 
for lettering. The smallest lettering on 8 X 10 inch figures should 
be no less than '/, inch high. Cross-section paper (white with 
black lines) may be used, but it should not have more than 4 lines 
per inch. If finer ruled paper is used, the major division lines 
should be drawn in with black ink, omitting the finer divisions. 
In the case of finely ruled paper, only blue-lined paper can be 
accepted. Tracing paper and blueprints are not acceptable. 
Lettering and all markings must be large enough to be readable 
after reduction to single-column width (35/:,in.). Mail rolled or 
flat; never fold. Drawings that cannot be reproduced (including 
pencil drawings) will be returned to the author for redrawing, 
thus delaying publication of the paper. Photographs should be 
distinct and show clear black and white contrasts. They must 
be on glossy white paper. Avoid round and oval photographs. 


CAPTIONS AND LEGENDS: Legends or captions must accompany 
each drawing or photograph submitted. If written on the draw- 
ing or photograph, they should be placed below and well outside 
the part to be reproduced. Each table should have a caption 
such as Table 1, Table 2, Table 3, etc. Captions should be com- 
plete in themselves so as to make the data intelligible to the 
reader without reference to the text. A duplicate list of captions 
for figures should be included as the last page of the manuscript. 
Use “Fig. 1’’ (not Figure 1), ‘‘Figs. 3 and 4,” etc., in both the text 
and the numbering of illustrations. In the text, ‘‘Eq. (1)”’ or 
“Eqs. (1) and is used, not “Equation In captions, 
legends, and in table headings, write all words in full; do not 
abbreviate, except for “Fig.” and ‘“‘Eq.” 


MATHEMATICAL WorRK: Formulas may be typewritten or 
carefully written in pen and ink, the writing to be large enough so 
that it can be marked for the printer. Considerable space for mark- 
ing should be allowed above and below all equations. All compli- 
cated equations should be repeated on separate sheets with plenty 
of space left for marking. The solidus should be used for sim- 
ple fractions appearing within the text. Make all expressions 
clear to the typesetter. Greek letters used in formulas should 
be clearly designated by name in the margin of the manu- 
script. The difference between capital and lower-case letters 
should be clearly distinguished and care taken to avoid confu- 
sion between zero (0) and the letter (0), between the numeral 
(one), and the letter (ell) and the prime (’), between alpha and 
a, kappa and k, u and mu, v and nu, n and eta. All sub- 
scripts and exponents should be clearly distinguished. Avoid 
complicated exponents and subscripts. Dots and bars over 
letters or mathematical expressions should also be avoided. 
When it is necessary to repeat a complicated expression, it should 
be represented by some convenient symbol. 


SYMBOLS AND ABBREVIATIONS: The symbols recommended in 
the American Standard Association ‘Letter Symbols for Aero- 
nautical Sciences,’”” ASA Y10.7—1954, should be used wherever 
practicable. All symbols should be clearly written and carefully 
checked. Standard abbreviations should be used, and it should 
be noted that most abbreviations are lower case, such as m.p.h., 
b.m.e.p., ib.p., b.hp., hp., etc. 


= é 
4 
| 
juz 
| 
¥ 
| 
=< 


CorRPORATE MEMBERS 


OF THE 


INSTITUTE OF THE AERONAUTICAL SCIENCES 


AC SPARK PLUG DIVISION, GENERAL MOTORS COR- 
PORATION 


ACADEMY OF AERONAUTICS 

ADAMS RITE MANUFACTURING COMPANY 
AEROJET-GENERAL CORPORATION 
AEROLAB DEVELOPMENT COMPANY 
AERONCA MANUFACTURING CORPORATION 
AERONUTRONIC SYSTEMS, INC. 

AEROQUIP CORPORATION 

AGAWAM AIRCRAFT PRODUCTS, INC. 
ALLIED RESEARCH ASSOCIATES, INC. 
ALLISON DIVISION, GENERAL MOTORS CORPORATION 
ALUMINUM COMPANY OF AMERICA 
AMERICAN AIRLINES, INC. 

AMERICAN BOSCH ARMA CORPORATION 


AMERICAN STEEL & WIRE DIVISION, UNITED STATES 
STEEL CORPORATION 


AMOCO CHEMICALS CORPORATION 
AMPHENOL-BORG ELECTRONICS CORPORATION 
ANDERSON, GREENWOOD & CO. 

ARDE ASSOCIATES 

ASSOCIATED AVIATION UNDERWRITERS 

AVCO RESEARCH LABORATORY 

AVIEN, INC. 

AVRO AIRCRAFT LIMITED 

BAKER STEEL & TUBE COMPANY 


BEECH AIRCRAFT CORPORATION 
BELL AIRCRAFT CORPORATION 


BENDIX AVIATION CORPORATION 

BOEING AIRPLANE COMPANY 

BOOZ, ALLEN & HAMILTON 

THE BRISTOL AEROPLANE COMPANY (U.S.A.) INC. 


BULOVA RESEARCH AND DEVELOPMENT LABORA- 
TORIES, INC. 


CANADAIR, LTD. 

THE CESSNA AIRCRAFT COMPANY 

CHANCE VOUGHT AIRCRAFT, INCORPORATED 
THE CHASE MANHATTAN BANK 

CHICAGO AERIAL INDUSTRIES, INC. 

THE CLEVELAND PNEUMATIC INDUSTRIES, INC. 
CONSOLIDATED CONTROLS CORPORATION 
CONTINENTAL MOTORS CORPORATION 
CORNELL AERONAUTICAL LABORATORY, INC. 
CURTISS-WRIGHT CORPORATION 

DANIEL, MANN, JOHNSON, & MENDENHALL 
THE DECKER CORPORATION 

DEL MAR ENGINEERING LABORATORIES 


DELCO-REMY DIVISION, GENERAL MOTORS COR- 
PORATION 


DOAK AIRCRAFT COMPANY, INC. 
DOUGLAS AIRCRAFT COMPANY, INC. 
DZUS FASTENER COMPANY, INC. 
EASTERN AIR LINES, INC. 


EATON MANUFACTURING COMPANY 

EDO CORPORATION 

ELASTIC STOP NUT CORPORATION OF AMERICA 

ELECTROL INCORPORATED 

ENGINEERING SUPERVISION COMPANY 

ESSO STANDARD OIL COMPANY 

FAIRCHILD CAMERA AND INSTRUMENT CORPORATION 

FAIRCHILD ENGINE AND AIRPLANE CORPORATION 

THE FIRESTONE TIRE & RUBBER COMPANY 

THE GARRETT CORPORATION 

GENERAL DYNAMICS CORPORATION 

GENERAL ELECTRIC COMPANY 

GENERAL PRECISION EQUIPMENT CORPORATION 

G. M. GIANNINI & CO., INC. 

GIANNINI PLASMADYNE CORPORATION 

THE B. F, GOODRICH COMPANY 

GOODYEAR AIRCRAFT CORPORATION 

GRUMMAN AIRCRAFT ENGINEERING CORPORATION 

HYDRO-AIRE, INC. 

INSURANCE COMPANY OF NORTH AMERICA COM- 
PANIES 


INTERNATIONAL BUSINESS MACHINES CORPORATION 
THE INTERNATIONAL NICKEL COMPANY, INC. 


ITT ae DIVISION OF INTERNATIONAL 
NE AND TELEGRAPH CORPORATION 


JANITROL AIRCRAFT DIVISION, SURFACE COMBUS- 
TION CORPORATION 


JOHNS-MANVILLE SALES CORPORATION 
EARLE M. JORGENSEN CO. 

WALTER KIDDE & COMPANY, INC. 
KOLLSMAN INSTRUMENT CORPORATION 
LAVELLE AIRCRAFT CORPORATION 
LEAR, INCORPORATED 

C. W. LEMMERMAN, INC. 

LIBRASCOPE, INCORPORATED 

THE LIQUIDOMETER CORPORATION 
LOCKHEED AIRCRAFT CORPORATION 


LOEWY-HYDROPRESS DIVISION OF BALDWIN-LIMA- 
HAMILTON CORPORATION 


LONGINES-WITTNAUER WATCH COMPANY, INC. 
MARQUARDT CORPORATION 

THE MARTIN COMPANY 

MARTIN STEEL PRODUCTS CORPORATION 
McDONNELL AIRCRAFT CORPORATION 

MELETRON CORPORATION 

MID-CONTINENT MANUFACTURING, INC. 
MINNEAPOLIS-HONEY WELL REGULATOR COMPANY 
NATIONAL CREDIT OFFICE, INC. 

NORTH AMERICAN AVIATION, INC. 


NORTHROP CORPORATION 

NUCLEAR DEVELOPMENT CORPORATION OF AMERICAM 

PAN AMERICAN WORLD AIRWAYS, INC, : 

THE RALPH M. PARSONS COMPANY q 

PESCO PRODUCTS DIVISION, BORG-WARNER CORs 
PORATION 


PHILLIPS PETROLEUM COMPANY 
PIASECKI AIRCRAFT CORPORATION 


RADIO CORPORATION OF AMERICA 
DEFENSE ELECTRONIC PRODUCTS 


REPUBLIC AVIATION CORPORATION 

ROHR AIRCRAFT CORPORATION 

PAUL ROSENBERG ASSOCIATES 

RYAN AERONAUTICAL COMPANY 
SANDBERG-SERRELL CORPORATION 
SHAFER BEARING DIVISION, CHAIN BELT COMPANY 
SHELL OIL COMPANY 

SIMMONDS AEROCESSORIES, INC. 

SOCONY MOBIL OIL COMPANY, INC. 
SOLAR AIRCRAFT COMPANY 

SOUTHWEST PRODUCTS CO, 

SPACE TECHNOLOGY LABORATORIES, INC. 
R. DIXON SPEAS ASSOCIATES 


SPERRY GYROSCOPE COMPANY DIVISION OF SPERRY¥ 
RAND CORPORATION 


STANDARD OIL COMPANY OF CALIFORNIA 
STANDARD OIL COMPANY (INDIANA) 


STANDARD-THOMSON CORPORATION 
CLIFFORD MANUFACTURING DIVISION 


STANLEY AVIATION CORPORATION 3 
THIOKOL CHEMICAL CORPORATION 
THOMPSON RAMO WOOLDRIDGE INC 
TRANS WORLD AIRLINES, INC. 
TURBO PRODUCTS, INC. 

UNION CARBIDE CORPORATION a 
UNITED AIR LINES, INC. 

UNITED AIRCRAFT CORPORATION 

UNITED STATES AVIATION UNDERWRITERS, INC 
VARD, INC. 

VERTOL AIRCRAFT CORPORATION 

VITRO CORPORATION OF AMERICA 

WESTERN GEAR CORPORATION 
WESTINGHOUSE ELECTRIC CORPORATION 
WESTON INSTRUMENTS DIVISION OF DAYSTROM& 


WYMAN-GORDON COMPANY 3 
YOUNG RADIATOR COMPANY q 


; 

it g 

a 

re 

i 

‘ 
‘ane 

q 

xs 

te 2 

fe 

io 


